From aec6e26c5fba23e8ddb43108aae7ff29710246a9 Mon Sep 17 00:00:00 2001 From: Dave Kleinschmidt Date: Mon, 21 Aug 2017 13:56:00 -0400 Subject: [PATCH 1/8] add a NormalInverseChisq prior for normal with unknown mean/variance --- src/ConjugatePriors.jl | 1 + src/normal.jl | 13 +++++++++ src/normalinversechisq.jl | 58 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 72 insertions(+) create mode 100644 src/normalinversechisq.jl diff --git a/src/ConjugatePriors.jl b/src/ConjugatePriors.jl index 0c750e7..5532efa 100644 --- a/src/ConjugatePriors.jl +++ b/src/ConjugatePriors.jl @@ -43,6 +43,7 @@ include("gamma_exp.jl") include("normalgamma.jl") include("normalinversegamma.jl") +include("normalinversechisq.jl") include("normalwishart.jl") include("normalinversewishart.jl") include("normal.jl") diff --git a/src/normal.jl b/src/normal.jl index be26ea6..5a1781d 100644 --- a/src/normal.jl +++ b/src/normal.jl @@ -111,6 +111,19 @@ end complete(G::Type{Normal}, pri::NormalInverseGamma, t::NTuple{2,Float64}) = Normal(t[1], sqrt(t[2])) +#### NormalInverseChisq on (μ, σ^2) + +function posterior_canon(prior::NormalInverseChisq, ss::NormalStats) + κn = prior.κ + ss.tw + νn = prior.ν + ss.tw + μn = (prior.κ*prior.μ + ss.s) / κn + σ2n = (prior.ν*prior.σ2 + ss.s2 + ss.tw*prior.κ/(ss.tw + prior.κ)*abs2(prior.μ - ss.m)) / νn + + NormalInverseChisq(μn, σ2n, κ, ν) +end + +complete(G::Type{Normal}, pri::NormalInverseChisq, t::NTuple{2,Float64}) = Normal(t[1], sqrt(t[2])) + #### NormalGamma on (μ, σ^(-2)) diff --git a/src/normalinversechisq.jl b/src/normalinversechisq.jl new file mode 100644 index 0000000..8fcc963 --- /dev/null +++ b/src/normalinversechisq.jl @@ -0,0 +1,58 @@ +# Based on Murphy (2007) Conjugate Bayesian Analysis of the Gaussian +# distribution. Equivalent to the Normal-Inverse Gamma under the follow +# re-parametrization: +# +# m0 = μ +# v0 = 1/κ +# shape = ν/2 +# scale = νσ2/2 +# +# The parameters have a natural interpretation when used as a prior for a Normal +# distribution with unknown mean and variance: μ and σ2 are the expected mean +# and variance, while κ and ν are the respective degrees of confidence +# (expressed in "pseudocounts"). + +struct NormalInverseChisq{T<:Real} <: ContinuousUnivariateDistribution + μ::T + σ2::T + κ::T + ν::T + + function NormalInverseChisq{T}(μ::T, σ2::T, κ::T, ν::T) where T<:Real + ν > zero(ν) && + κ > zero(κ) && + σ2 > zero(σ2) || + error("Variance and confidence (κ and ν) must all be positive") + new{T}(μ, σ2, κ, ν) + end +end + +function NormalInverseChisq(μ::Real, σ2::Real, κ::Real, ν::Real) + T = promote_type(typeof(μ), typeof(σ2), typeof(κ), typeof(ν)) + NormalInverseChisq{T}(T(μ), T(σ2), T(κ), T(ν)) +end + +Base.convert(::Type{NormalInverseGamma}, d::NormalInverseChisq) = + NormalInverseGamma(d.μ, 1/d.κ, d.ν/2, d.ν*d.σ2/2) + +Base.convert(::Type{NormalInverseChisq}, d::NormalInverseGamma) = + NormalInverseChisq(d.mu, d.scale/d.shape, 1/d.v0, d.shape*2) + +insupport(::Type{NormalInverseChisq}, μ::T, σ2::T) where T<:Real = + isfinite(μ) && zero(σ2) <= σ2 < Inf + +# function pdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real +# Zinv = sqrt(d.κ / 2pi) / gamma(d.ν*0.5) * (d.ν * d.σ2 / 2)^(d.ν*0.5) +# Zinv * σ2^(-(d.ν+3)*0.5) * exp( (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2)) +# end + +# function logpdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real +# logZinv = (log(d.κ) - log(2pi))*0.5 - lgamma(d.ν*0.5) + (log(d.ν) + log(d.σ2) - log(2)) * (d.ν/2) +# logZinv + log(σ2)*(-(d.ν+3)*0.5) + (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2) +# end + +pdf(d::NormalInverseChisq, μ::T, σ2::T) = pdf(NormalInverseGamma(d), μ, σ2) +logpdf(d::NormalInverseChisq, μ::T, σ2::T) = logpdf(NormalInverseGamma(d), μ, σ2) +mean(d::NormalInverseChisq) = mean(NormalInverseGamma(d)) +mode(d::NormalInverseChisq) = mode(NormalInverseGamma(d)) +rand(d::NormalInverseChisq) = mean(NormalInverseGamma(d)) From ebfc374707d4106a91676f26a83a293b6ca2764a Mon Sep 17 00:00:00 2001 From: Dave Kleinschmidt Date: Mon, 21 Aug 2017 14:20:21 -0400 Subject: [PATCH 2/8] fix bugs and update tests for normals to use testset --- src/normal.jl | 2 +- src/normalinversechisq.jl | 6 +- test/conjugates_normal.jl | 255 ++++++++++++++++++++++---------------- 3 files changed, 151 insertions(+), 112 deletions(-) diff --git a/src/normal.jl b/src/normal.jl index 5a1781d..11bc432 100644 --- a/src/normal.jl +++ b/src/normal.jl @@ -119,7 +119,7 @@ function posterior_canon(prior::NormalInverseChisq, ss::NormalStats) μn = (prior.κ*prior.μ + ss.s) / κn σ2n = (prior.ν*prior.σ2 + ss.s2 + ss.tw*prior.κ/(ss.tw + prior.κ)*abs2(prior.μ - ss.m)) / νn - NormalInverseChisq(μn, σ2n, κ, ν) + NormalInverseChisq(μn, σ2n, κn, νn) end complete(G::Type{Normal}, pri::NormalInverseChisq, t::NTuple{2,Float64}) = Normal(t[1], sqrt(t[2])) diff --git a/src/normalinversechisq.jl b/src/normalinversechisq.jl index 8fcc963..82e0e02 100644 --- a/src/normalinversechisq.jl +++ b/src/normalinversechisq.jl @@ -51,8 +51,8 @@ insupport(::Type{NormalInverseChisq}, μ::T, σ2::T) where T<:Real = # logZinv + log(σ2)*(-(d.ν+3)*0.5) + (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2) # end -pdf(d::NormalInverseChisq, μ::T, σ2::T) = pdf(NormalInverseGamma(d), μ, σ2) -logpdf(d::NormalInverseChisq, μ::T, σ2::T) = logpdf(NormalInverseGamma(d), μ, σ2) +pdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real = pdf(NormalInverseGamma(d), μ, σ2) +logpdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real = logpdf(NormalInverseGamma(d), μ, σ2) mean(d::NormalInverseChisq) = mean(NormalInverseGamma(d)) mode(d::NormalInverseChisq) = mode(NormalInverseGamma(d)) -rand(d::NormalInverseChisq) = mean(NormalInverseGamma(d)) +rand(d::NormalInverseChisq) = rand(NormalInverseGamma(d)) diff --git a/test/conjugates_normal.jl b/test/conjugates_normal.jl index e2af720..9e48ed2 100644 --- a/test/conjugates_normal.jl +++ b/test/conjugates_normal.jl @@ -1,158 +1,197 @@ -# Conjugates for normal distribution +# using Base.Test using Distributions using ConjugatePriors -import ConjugatePriors: NormalGamma, NormalInverseGamma -import ConjugatePriors: posterior, posterior_rand, posterior_mode, posterior_randmodel, fit_map +using ConjugatePriors: NormalGamma, NormalInverseGamma, NormalInverseChisq +using ConjugatePriors: posterior, posterior_rand, posterior_mode, posterior_randmodel, fit_map n = 100 w = rand(100) -# Νormal - Νormal (known sigma) +@testset "Conjugates for normal distribution" begin -pri = Normal(1.0, 5.0) + @testset "Νormal - Νormal (known sigma)" begin -x = rand(Normal(2.0, 3.0), n) -p = posterior((pri, 3.0), Normal, x) -@test isa(p, Normal) -@test mean(p) ≈ (mean(pri) / var(pri) + sum(x) / 9.0) / (1.0 / var(pri) + n / 9.0) -@test var(p) ≈ inv(1.0 / var(pri) + n / 9.0) + pri = Normal(1.0, 5.0) -r = posterior_mode((pri, 3.0), Normal, x) -@test r ≈ mode(p) + x = rand(Normal(2.0, 3.0), n) + p = posterior((pri, 3.0), Normal, x) + @test isa(p, Normal) + @test mean(p) ≈ (mean(pri) / var(pri) + sum(x) / 9.0) / (1.0 / var(pri) + n / 9.0) + @test var(p) ≈ inv(1.0 / var(pri) + n / 9.0) -f = fit_map((pri, 3.0), Normal, x) -@test isa(f, Normal) -@test f.μ == r -@test f.σ == 3.0 + r = posterior_mode((pri, 3.0), Normal, x) + @test r ≈ mode(p) -p = posterior((pri, 3.0), Normal, x, w) -@test isa(p, Normal) -@test mean(p) ≈ (mean(pri) / var(pri) + dot(x, w) / 9.0) / (1.0 / var(pri) + sum(w) / 9.0) -@test var(p) ≈ inv(1.0 / var(pri) + sum(w) / 9.0) + f = fit_map((pri, 3.0), Normal, x) + @test isa(f, Normal) + @test f.μ == r + @test f.σ == 3.0 -r = posterior_mode((pri, 3.0), Normal, x, w) -@test r ≈ mode(p) + p = posterior((pri, 3.0), Normal, x, w) + @test isa(p, Normal) + @test mean(p) ≈ (mean(pri) / var(pri) + dot(x, w) / 9.0) / (1.0 / var(pri) + sum(w) / 9.0) + @test var(p) ≈ inv(1.0 / var(pri) + sum(w) / 9.0) -f = fit_map((pri, 3.0), Normal, x, w) -@test isa(f, Normal) -@test f.μ == r -@test f.σ == 3.0 + r = posterior_mode((pri, 3.0), Normal, x, w) + @test r ≈ mode(p) + f = fit_map((pri, 3.0), Normal, x, w) + @test isa(f, Normal) + @test f.μ == r + @test f.σ == 3.0 -# ΙnverseGamma - Νormal (known mu) + end -pri = InverseGamma(1.5, 0.5) + @testset "ΙnverseGamma - Νormal (known mu)" begin -x = rand(Normal(2.0, 3.0), n) -p = posterior((2.0, pri), Normal, x) -@test isa(p, InverseGamma) -@test shape(p) ≈ shape(pri) + n / 2 -@test scale(p) ≈ scale(pri) + sum(abs2.(x .- 2.0)) / 2 + pri = InverseGamma(1.5, 0.5) -r = posterior_mode((2.0, pri), Normal, x) -@test r ≈ mode(p) + x = rand(Normal(2.0, 3.0), n) + p = posterior((2.0, pri), Normal, x) + @test isa(p, InverseGamma) + @test shape(p) ≈ shape(pri) + n / 2 + @test scale(p) ≈ scale(pri) + sum(abs2.(x .- 2.0)) / 2 -f = fit_map((2.0, pri), Normal, x) -@test isa(f, Normal) -@test f.μ == 2.0 -@test abs2(f.σ) ≈ r + r = posterior_mode((2.0, pri), Normal, x) + @test r ≈ mode(p) -p = posterior((2.0, pri), Normal, x, w) -@test isa(p, InverseGamma) -@test shape(p) ≈ shape(pri) + sum(w) / 2 -@test scale(p) ≈ scale(pri) + dot(w, abs2.(x .- 2.0)) / 2 + f = fit_map((2.0, pri), Normal, x) + @test isa(f, Normal) + @test f.μ == 2.0 + @test abs2(f.σ) ≈ r -r = posterior_mode((2.0, pri), Normal, x, w) -@test r ≈ mode(p) + p = posterior((2.0, pri), Normal, x, w) + @test isa(p, InverseGamma) + @test shape(p) ≈ shape(pri) + sum(w) / 2 + @test scale(p) ≈ scale(pri) + dot(w, abs2.(x .- 2.0)) / 2 -f = fit_map((2.0, pri), Normal, x, w) -@test isa(f, Normal) -@test f.μ == 2.0 -@test abs2(f.σ) ≈ r + r = posterior_mode((2.0, pri), Normal, x, w) + @test r ≈ mode(p) + f = fit_map((2.0, pri), Normal, x, w) + @test isa(f, Normal) + @test f.μ == 2.0 + @test abs2(f.σ) ≈ r -# Gamma - Νormal (known mu) + end -pri = Gamma(1.5, 2.0) + @testset "Gamma - Νormal (known mu)" begin -x = rand(Normal(2.0, 3.0), n) -p = posterior((2.0, pri), Normal, x) -@test isa(p, Gamma) -@test shape(p) ≈ shape(pri) + n / 2 -@test scale(p) ≈ scale(pri) + sum(abs2.(x .- 2.0)) / 2 + pri = Gamma(1.5, 2.0) -r = posterior_mode((2.0, pri), Normal, x) -@test r ≈ mode(p) + x = rand(Normal(2.0, 3.0), n) + p = posterior((2.0, pri), Normal, x) + @test isa(p, Gamma) + @test shape(p) ≈ shape(pri) + n / 2 + @test scale(p) ≈ scale(pri) + sum(abs2.(x .- 2.0)) / 2 -f = fit_map((2.0, pri), Normal, x) -@test isa(f, Normal) -@test f.μ == 2.0 -@test abs2(f.σ) ≈ inv(r) + r = posterior_mode((2.0, pri), Normal, x) + @test r ≈ mode(p) -p = posterior((2.0, pri), Normal, x, w) -@test isa(p, Gamma) -@test shape(p) ≈ shape(pri) + sum(w) / 2 -@test scale(p) ≈ scale(pri) + dot(w, abs2.(x .- 2.0)) / 2 + f = fit_map((2.0, pri), Normal, x) + @test isa(f, Normal) + @test f.μ == 2.0 + @test abs2(f.σ) ≈ inv(r) -r = posterior_mode((2.0, pri), Normal, x, w) -@test r ≈ mode(p) + p = posterior((2.0, pri), Normal, x, w) + @test isa(p, Gamma) + @test shape(p) ≈ shape(pri) + sum(w) / 2 + @test scale(p) ≈ scale(pri) + dot(w, abs2.(x .- 2.0)) / 2 -f = fit_map((2.0, pri), Normal, x, w) -@test isa(f, Normal) -@test f.μ == 2.0 -@test abs2(f.σ) ≈ inv(r) + r = posterior_mode((2.0, pri), Normal, x, w) + @test r ≈ mode(p) + f = fit_map((2.0, pri), Normal, x, w) + @test isa(f, Normal) + @test f.μ == 2.0 + @test abs2(f.σ) ≈ inv(r) -# NormalInverseGamma - Normal + end -mu_true = 2. -sig2_true = 3. -x = rand(Normal(mu_true, sig2_true), n) + @testset "NormalInverseGamma - Normal" begin -mu0 = 2. -v0 = 3. -shape0 = 5. -scale0 = 2. -pri = NormalInverseGamma(mu0, v0, shape0, scale0) + mu_true = 2. + sig2_true = 3. + x = rand(Normal(mu_true, sig2_true), n) -post = posterior(pri, Normal, x) -@test isa(post, NormalInverseGamma) + mu0 = 2. + v0 = 3. + shape0 = 5. + scale0 = 2. + pri = NormalInverseGamma(mu0, v0, shape0, scale0) -@test post.mu ≈ (mu0/v0 + n*mean(x))/(1./v0 + n) -@test post.v0 ≈ 1./(1./v0 + n) -@test post.shape ≈ shape0 + 0.5*n -@test post.scale ≈ scale0 + 0.5*(n-1)*var(x) + n./v0./(n + 1./v0)*0.5*(mean(x)-mu0).^2 + post = posterior(pri, Normal, x) + @test isa(post, NormalInverseGamma) -ps = posterior_randmodel(pri, Normal, x) + @test post.mu ≈ (mu0/v0 + n*mean(x))/(1./v0 + n) + @test post.v0 ≈ 1./(1./v0 + n) + @test post.shape ≈ shape0 + 0.5*n + @test post.scale ≈ scale0 + 0.5*(n-1)*var(x) + n./v0./(n + 1./v0)*0.5*(mean(x)-mu0).^2 -@test isa(ps, Normal) -@test insupport(ps,ps.μ) && ps.σ > zero(ps.σ) + ps = posterior_randmodel(pri, Normal, x) + @test isa(ps, Normal) + @test insupport(ps,ps.μ) && ps.σ > zero(ps.σ) -# NormalGamma - Normal + end -mu_true = 2. -tau2_true = 3. -x = rand(Normal(mu_true, 1./tau2_true), n) + @testset "NormalInverseChisq - Normal" begin -mu0 = 2. -nu0 = 3. -shape0 = 5. -rate0 = 2. -pri = NormalGamma(mu0, nu0, shape0, rate0) + mu_true = 2. + sig2_true = 3. + x = rand(Normal(mu_true, sig2_true), n) -post = posterior(pri, Normal, x) -@test isa(post, NormalGamma) + μ0 = 2.0 + σ20 = 2.0/5.0 + κ0 = 1.0/3.0 + ν0 = 10.0 -@test post.mu ≈ (nu0*mu0 + n*mean(x))./(nu0 + n) -@test post.nu ≈ nu0 + n -@test post.shape ≈ shape0 + 0.5*n -@test post.rate ≈ rate0 + 0.5*(n-1)*var(x) + n*nu0/(n + nu0)*0.5*(mean(x)-mu0).^2 + pri = NormalInverseChisq(μ0, σ20, κ0, ν0) + pri2 = NormalInverseGamma(pri) -ps = posterior_randmodel(pri, Normal, x) + @test NormalInverseChisq(pri2) == pri -@test isa(ps, Normal) -@test insupport(ps, ps.μ) && ps.σ > zero(ps.σ) + @test mode(pri2) == mode(pri) + @test mean(pri2) == mean(pri) + @test pdf(pri, mu_true, sig2_true) == pdf(pri2, mu_true, sig2_true) + + @test (srand(1); rand(pri)) == (srand(1); rand(pri2)) + + # check that updating is consistent between NIχ2 and NIG + post = posterior(pri, Normal, x) + post2 = posterior(pri2, Normal, x) + @test isa(post, NormalInverseChisq) + @test NormalInverseChisq(post2) == post + + end + + @testset "NormalGamma - Normal" begin + + mu_true = 2. + tau2_true = 3. + x = rand(Normal(mu_true, 1./tau2_true), n) + + mu0 = 2. + nu0 = 3. + shape0 = 5. + rate0 = 2. + pri = NormalGamma(mu0, nu0, shape0, rate0) + + post = posterior(pri, Normal, x) + @test isa(post, NormalGamma) + + @test post.mu ≈ (nu0*mu0 + n*mean(x))./(nu0 + n) + @test post.nu ≈ nu0 + n + @test post.shape ≈ shape0 + 0.5*n + @test post.rate ≈ rate0 + 0.5*(n-1)*var(x) + n*nu0/(n + nu0)*0.5*(mean(x)-mu0).^2 + + ps = posterior_randmodel(pri, Normal, x) + + @test isa(ps, Normal) + @test insupport(ps, ps.μ) && ps.σ > zero(ps.σ) + end + +end From 607af6db758d95426fc29dc44fbfc9268d058506 Mon Sep 17 00:00:00 2001 From: Dave Kleinschmidt Date: Mon, 4 Sep 2017 15:01:14 -0400 Subject: [PATCH 3/8] params method and no-params constructor for NormalInverseChisq --- src/ConjugatePriors.jl | 5 ++++- src/normalinversechisq.jl | 8 ++++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/ConjugatePriors.jl b/src/ConjugatePriors.jl index 5532efa..7f704e1 100644 --- a/src/ConjugatePriors.jl +++ b/src/ConjugatePriors.jl @@ -8,6 +8,8 @@ using Distributions import Base: mean import Base.LinAlg: Cholesky + + import Distributions: BernoulliStats, BinomData, @@ -34,7 +36,8 @@ import Distributions: mode, rand, pdf, - logpdf + logpdf, + params include("fallbacks.jl") include("beta_binom.jl") diff --git a/src/normalinversechisq.jl b/src/normalinversechisq.jl index 82e0e02..5454ec6 100644 --- a/src/normalinversechisq.jl +++ b/src/normalinversechisq.jl @@ -19,14 +19,16 @@ struct NormalInverseChisq{T<:Real} <: ContinuousUnivariateDistribution ν::T function NormalInverseChisq{T}(μ::T, σ2::T, κ::T, ν::T) where T<:Real - ν > zero(ν) && - κ > zero(κ) && + ν >= zero(ν) && + κ >= zero(κ) && σ2 > zero(σ2) || error("Variance and confidence (κ and ν) must all be positive") new{T}(μ, σ2, κ, ν) end end +NormalInverseChisq() = NormalInverseChisq{Float64}(0.0, 1.0, 0.0, 0.0) + function NormalInverseChisq(μ::Real, σ2::Real, κ::Real, ν::Real) T = promote_type(typeof(μ), typeof(σ2), typeof(κ), typeof(ν)) NormalInverseChisq{T}(T(μ), T(σ2), T(κ), T(ν)) @@ -41,6 +43,8 @@ Base.convert(::Type{NormalInverseChisq}, d::NormalInverseGamma) = insupport(::Type{NormalInverseChisq}, μ::T, σ2::T) where T<:Real = isfinite(μ) && zero(σ2) <= σ2 < Inf +params(d::NormalInverseChisq) = d.μ, d.σ2, d.κ, d.ν + # function pdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real # Zinv = sqrt(d.κ / 2pi) / gamma(d.ν*0.5) * (d.ν * d.σ2 / 2)^(d.ν*0.5) # Zinv * σ2^(-(d.ν+3)*0.5) * exp( (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2)) From 1364ba836ffc55c16c4e0857924712bddb26b881 Mon Sep 17 00:00:00 2001 From: Dave Kleinschmidt Date: Tue, 5 Dec 2017 20:54:45 -0500 Subject: [PATCH 4/8] docstring for nix2 and fix conditional --- src/normalinversechisq.jl | 45 ++++++++++++++++++++++++--------------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/src/normalinversechisq.jl b/src/normalinversechisq.jl index 5454ec6..2a3948c 100644 --- a/src/normalinversechisq.jl +++ b/src/normalinversechisq.jl @@ -1,17 +1,29 @@ -# Based on Murphy (2007) Conjugate Bayesian Analysis of the Gaussian -# distribution. Equivalent to the Normal-Inverse Gamma under the follow -# re-parametrization: -# -# m0 = μ -# v0 = 1/κ -# shape = ν/2 -# scale = νσ2/2 -# -# The parameters have a natural interpretation when used as a prior for a Normal -# distribution with unknown mean and variance: μ and σ2 are the expected mean -# and variance, while κ and ν are the respective degrees of confidence -# (expressed in "pseudocounts"). +""" + NormalInverseChisq(μ, σ2, κ, ν) +A Normal-χ^-2 distribution is a conjugate prior for a Normal distribution with +unknown mean and variance. It has parameters: + +* μ: expected mean +* σ2 > 0: expected variance +* κ ≥ 0: mean confidence +* ν ≥ 0: variance confidence + +The parameters have a natural interpretation when used as a prior for a Normal +distribution with unknown mean and variance: μ and σ2 are the expected mean and +variance, while κ and ν are the respective degrees of confidence (expressed in +"pseudocounts"). When interpretable parameters are important, this makes it a +slightly more convient parametrization of the conjugate prior. + +Equivalent to a Normal-Inverse Gamma distribution with parameters: + +* m0 = μ +* v0 = 1/κ +* shape = ν/2 +* scale = νσ2/2 + +Based on Murphy "Conjugate Bayesian analysis of the Gaussian distribution". +""" struct NormalInverseChisq{T<:Real} <: ContinuousUnivariateDistribution μ::T σ2::T @@ -19,10 +31,9 @@ struct NormalInverseChisq{T<:Real} <: ContinuousUnivariateDistribution ν::T function NormalInverseChisq{T}(μ::T, σ2::T, κ::T, ν::T) where T<:Real - ν >= zero(ν) && - κ >= zero(κ) && - σ2 > zero(σ2) || - error("Variance and confidence (κ and ν) must all be positive") + if ν < 0 || κ < 0 || σ2 ≤ 0 + throw(ArgumentError("Variance and confidence (κ and ν) must all be positive")) + end new{T}(μ, σ2, κ, ν) end end From abb86686284c4da1be1a1ea6d6c43c063c20b13b Mon Sep 17 00:00:00 2001 From: Dave Kleinschmidt Date: Tue, 5 Dec 2017 21:49:56 -0500 Subject: [PATCH 5/8] use native pdf, logpdf, mean, mode functions --- src/normalinversechisq.jl | 36 ++++++++++++++++++++++-------------- test/conjugates_normal.jl | 6 ++++++ 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/src/normalinversechisq.jl b/src/normalinversechisq.jl index 2a3948c..8342a86 100644 --- a/src/normalinversechisq.jl +++ b/src/normalinversechisq.jl @@ -56,18 +56,26 @@ insupport(::Type{NormalInverseChisq}, μ::T, σ2::T) where T<:Real = params(d::NormalInverseChisq) = d.μ, d.σ2, d.κ, d.ν -# function pdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real -# Zinv = sqrt(d.κ / 2pi) / gamma(d.ν*0.5) * (d.ν * d.σ2 / 2)^(d.ν*0.5) -# Zinv * σ2^(-(d.ν+3)*0.5) * exp( (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2)) -# end - -# function logpdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real -# logZinv = (log(d.κ) - log(2pi))*0.5 - lgamma(d.ν*0.5) + (log(d.ν) + log(d.σ2) - log(2)) * (d.ν/2) -# logZinv + log(σ2)*(-(d.ν+3)*0.5) + (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2) -# end - -pdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real = pdf(NormalInverseGamma(d), μ, σ2) -logpdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real = logpdf(NormalInverseGamma(d), μ, σ2) -mean(d::NormalInverseChisq) = mean(NormalInverseGamma(d)) -mode(d::NormalInverseChisq) = mode(NormalInverseGamma(d)) +function pdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real + Zinv = sqrt(d.κ / 2pi) / gamma(d.ν*0.5) * (d.ν * d.σ2 / 2)^(d.ν*0.5) + Zinv * σ2^(-(d.ν+3)*0.5) * exp( (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2)) +end + +function logpdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real + logZinv = (log(d.κ) - log(2pi))*0.5 - lgamma(d.ν*0.5) + (log(d.ν) + log(d.σ2) - log(2)) * (d.ν/2) + logZinv + log(σ2)*(-(d.ν+3)*0.5) + (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2) +end + +function mean(d::NormalInverseChisq) + μ = d.μ + σ2 = d.ν/(d.ν-2)*d.σ2 + return μ, σ2 +end + +function mode(d::NormalInverseChisq) + μ = d.μ + σ2 = d.ν*d.σ2/(d.ν - 1) + return μ, σ2 +end + rand(d::NormalInverseChisq) = rand(NormalInverseGamma(d)) diff --git a/test/conjugates_normal.jl b/test/conjugates_normal.jl index 9e48ed2..43cb02b 100644 --- a/test/conjugates_normal.jl +++ b/test/conjugates_normal.jl @@ -166,6 +166,12 @@ w = rand(100) @test isa(post, NormalInverseChisq) @test NormalInverseChisq(post2) == post + for _ in 1:10 + x = rand(post) + @test pdf(post, x...) ≈ pdf(post2, x...) + @test logpdf(post, x...) ≈ logpdf(post2, x...) + end + end @testset "NormalGamma - Normal" begin From 98fe3948b6b9428128108e59732004a96de6e829 Mon Sep 17 00:00:00 2001 From: Dave Kleinschmidt Date: Wed, 6 Dec 2017 09:37:25 -0500 Subject: [PATCH 6/8] fix docstring: typo and refer to NormalInverseGamma --- src/normalinversechisq.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/normalinversechisq.jl b/src/normalinversechisq.jl index 8342a86..4ff6e7e 100644 --- a/src/normalinversechisq.jl +++ b/src/normalinversechisq.jl @@ -13,9 +13,9 @@ The parameters have a natural interpretation when used as a prior for a Normal distribution with unknown mean and variance: μ and σ2 are the expected mean and variance, while κ and ν are the respective degrees of confidence (expressed in "pseudocounts"). When interpretable parameters are important, this makes it a -slightly more convient parametrization of the conjugate prior. +slightly more convenient parametrization of the conjugate prior. -Equivalent to a Normal-Inverse Gamma distribution with parameters: +Equivalent to a `NormalInverseGamma` distribution with parameters: * m0 = μ * v0 = 1/κ From d58befe650d27c6ca4e1c28ebe0ed568fd6f54f3 Mon Sep 17 00:00:00 2001 From: Dave Kleinschmidt Date: Wed, 6 Dec 2017 09:45:15 -0500 Subject: [PATCH 7/8] modes are not equal --- test/conjugates_normal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/conjugates_normal.jl b/test/conjugates_normal.jl index 43cb02b..6ed7a85 100644 --- a/test/conjugates_normal.jl +++ b/test/conjugates_normal.jl @@ -154,7 +154,7 @@ w = rand(100) @test NormalInverseChisq(pri2) == pri - @test mode(pri2) == mode(pri) + @test_broken mode(pri2) == mode(pri) @test mean(pri2) == mean(pri) @test pdf(pri, mu_true, sig2_true) == pdf(pri2, mu_true, sig2_true) From 77f0c6429326af31c7f19a0af5820469de1ac5bd Mon Sep 17 00:00:00 2001 From: Dave Kleinschmidt Date: Wed, 6 Dec 2017 09:46:10 -0500 Subject: [PATCH 8/8] remove unnecessary whitespace change --- src/ConjugatePriors.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ConjugatePriors.jl b/src/ConjugatePriors.jl index 7f704e1..53116c0 100644 --- a/src/ConjugatePriors.jl +++ b/src/ConjugatePriors.jl @@ -8,8 +8,6 @@ using Distributions import Base: mean import Base.LinAlg: Cholesky - - import Distributions: BernoulliStats, BinomData,