Merge pull request #18 from JuliaStats/dfk/normalinversechisq

Add normal-inverse Chi-squared parametrization

src/ConjugatePriors.jl

@@ -34,7 +34,8 @@ import Distributions:
     mode,
     rand,
     pdf,
-    logpdf
+    logpdf,
+    params
 
 include("fallbacks.jl")
 include("beta_binom.jl")
@@ -43,6 +44,7 @@ include("gamma_exp.jl")
 
 include("normalgamma.jl")
 include("normalinversegamma.jl")
+include("normalinversechisq.jl")
 include("normalwishart.jl")
 include("normalinversewishart.jl")
 include("normal.jl")

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, κn, νn)
+end
+
+complete(G::Type{Normal}, pri::NormalInverseChisq, t::NTuple{2,Float64}) = Normal(t[1], sqrt(t[2]))
+
 
 #### NormalGamma on (μ, σ^(-2))
 

src/normalinversechisq.jl

@@ -0,0 +1,81 @@
+"""
+    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 convenient parametrization of the conjugate prior.
+
+Equivalent to a `NormalInverseGamma` 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
+    κ::T
+    ν::T
+
+    function NormalInverseChisq{T}(μ::T, σ2::T, κ::T, ν::T) where T<:Real
+        if ν < 0 || κ < 0 || σ2 ≤ 0
+            throw(ArgumentError("Variance and confidence (κ and ν) must all be positive"))
+        end
+        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(ν))
+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
+
+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
+
+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))