Merge pull request #21 from JuliaStats/dfk/norminvgamma-mode

Use correct formula for normal-inverse chi2 mode and normal-inverse gamma mode
JuliaStats · Feb 13, 2018 · 8e3a6a7 · 8e3a6a7
2 parents 36fe174 + 0b1ec09
commit 8e3a6a7
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 3 deletions.
diff --git a/src/normalinversechisq.jl b/src/normalinversechisq.jl
@@ -74,7 +74,7 @@ end
 
 function mode(d::NormalInverseChisq)
     μ = d.μ
-    σ2 = d.ν*d.σ2/(d.ν - 1)
+    σ2 = d.ν*d.σ2/(d.ν + 3)
     return μ, σ2
 end
 

diff --git a/src/normalinversegamma.jl b/src/normalinversegamma.jl
@@ -43,7 +43,7 @@ end
 
 function mode(d::NormalInverseGamma)
     mu = d.mu
-    sig2 = d.scale / (d.shape + 1.0)
+    sig2 = d.scale / (d.shape + 1.5)
     return mu, sig2
 end
 

diff --git a/test/conjugates_normal.jl b/test/conjugates_normal.jl
@@ -133,6 +133,13 @@ w = rand(100)
         @test isa(ps, Normal)
         @test insupport(ps,ps.μ) && ps.σ > zero(ps.σ)
 
+        for shape in [0.1, 1., 10.]
+            d = NormalInverseGamma(mu0, v0, scale0, shape)
+            μ, σ2 = mode(d)
+            @test pdf(d, μ, σ2) > pdf(d, μ, σ2 + 0.001)
+            @test pdf(d, μ, σ2) > pdf(d, μ, σ2 - 0.001)
+        end
+
     end
 
     @testset "NormalInverseChisq - Normal" begin
@@ -151,7 +158,7 @@ w = rand(100)
 
         @test NormalInverseChisq(pri2) == pri
 
-        @test_broken mode(pri2) == mode(pri)
+        @test mode(pri2) == mode(pri)
         @test mean(pri2) == mean(pri)
         @test pdf(pri, mu_true, sig2_true) == pdf(pri2, mu_true, sig2_true)
 
@@ -169,6 +176,13 @@ w = rand(100)
             @test logpdf(post, x...) ≈ logpdf(post2, x...)
         end
 
+        for ν in [0.1, 1., 10.]
+            nix = NormalInverseChisq(0., 2., 3., ν)
+            μ, σ2 = mode(nix)
+            @test pdf(nix, μ, σ2) > pdf(nix, μ, σ2 + 0.001)
+            @test pdf(nix, μ, σ2) > pdf(nix, μ, σ2 - 0.001)
+        end
+
     end
 
     @testset "NormalGamma - Normal" begin