Commit a8d84af3 authored by Dahua Lin's avatar Dahua Lin

Change parameterization of InverseGamma (using the wikipedia way)

parent 1beadf17
......@@ -188,8 +188,8 @@ end
function posterior(prior::InverseGamma, ss::NormalKnownMuStats)
α1 = prior.shape + ss.tw / 2
β1 = rate(prior) + ss.s2 / 2
return InverseGamma(α1, 1.0 / β1)
β1 = prior.scale + ss.s2 / 2
return InverseGamma(α1, β1)
end
function posterior{T<:Real}(prior::(Float64, InverseGamma), ::Type{Normal}, x::Array{T})
......
##############################################################################
#
# InverseGamma distribution, cf Bayesian Theory (Barnardo & Smith) p 119
#
# Note that B&S parametrize in terms of shape/rate, but this is uses
# shape/scale so that it is consistent with the implementation of Gamma
#
##############################################################################
# Inverse Gamma distribution
immutable InverseGamma <: ContinuousUnivariateDistribution
shape::Float64 # location
......@@ -16,58 +9,63 @@ immutable InverseGamma <: ContinuousUnivariateDistribution
end
end
_inv(d::InverseGamma) = Gamma(d.shape, 1.0 / d.scale)
scale(d::InverseGamma) = d.scale
rate(d::InverseGamma) = 1.0 / d.scale
insupport(::InverseGamma, x::Real) = zero(x) <= x < Inf
insupport(::Type{InverseGamma}, x::Real) = zero(x) <= x < Inf
mean(d::InverseGamma) = d.shape > 1.0 ? (1.0 / d.scale) / (d.shape - 1.0) : Inf
var(d::InverseGamma) = d.shape > 2.0 ? (1.0 / d.scale)^2 / ((d.shape - 1.0)^2 * (d.shape - 2.0)) : Inf
mean(d::InverseGamma) = d.shape > 1.0 ? d.scale / (d.shape - 1.0) : Inf
var(d::InverseGamma) = d.shape > 2.0 ? abs2(d.scale) / (abs2(d.shape - 1.0) * (d.shape - 2.0)) : Inf
skewness(d::InverseGamma) = d.shape > 3.0 ? (4.0 * sqrt(d.shape - 2.0)) / (d.shape - 3.0) : NaN
kurtosis(d::InverseGamma) = d.shape > 4.0 ? (30.0 * d.shape - 66.0) / ((d.shape - 3.0) * (d.shape - 4.0)) : NaN
mode(d::InverseGamma) = (1.0 / d.scale) / (d.shape + 1.0)
mode(d::InverseGamma) = d.scale / (d.shape + 1.0)
modes(d::InverseGamma) = [mode(d)]
cdf(d::InverseGamma, x::Real) = ccdf(Gamma(d.shape, d.scale), 1.0 / x)
ccdf(d::InverseGamma, x::Real) = cdf(Gamma(d.shape, d.scale), 1.0 / x)
logcdf(d::InverseGamma, x::Real) = logccdf(Gamma(d.shape, d.scale), 1.0 / x)
logccdf(d::InverseGamma, x::Real) = logcdf(Gamma(d.shape, d.scale), 1.0 / x)
cdf(d::InverseGamma, x::Real) = ccdf(_inv(d), 1.0 / x)
ccdf(d::InverseGamma, x::Real) = cdf(_inv(d), 1.0 / x)
logcdf(d::InverseGamma, x::Real) = logccdf(_inv(d), 1.0 / x)
logccdf(d::InverseGamma, x::Real) = logcdf(_inv(d), 1.0 / x)
quantile(d::InverseGamma, p::Real) = 1.0 / cquantile(Gamma(d.shape, d.scale), p)
cquantile(d::InverseGamma, p::Real) = 1.0 / quantile(Gamma(d.shape, d.scale), p)
invlogcdf(d::InverseGamma, p::Real) = 1.0 / invlogccdf(Gamma(d.shape, d.scale), p)
invlogccdf(d::InverseGamma, p::Real) = 1.0 / invlogcdf(Gamma(d.shape, d.scale), p)
quantile(d::InverseGamma, p::Real) = 1.0 / cquantile(_inv(d), p)
cquantile(d::InverseGamma, p::Real) = 1.0 / quantile(_inv(d), p)
invlogcdf(d::InverseGamma, p::Real) = 1.0 / invlogccdf(_inv(d), p)
invlogccdf(d::InverseGamma, p::Real) = 1.0 / invlogcdf(_inv(d), p)
function entropy(d::InverseGamma)
a, b = d.shape, d.scale
a - log(b) + lgamma(a) - (1.0 + a) * digamma(a)
a = d.shape
b = d.scale
a + log(b) + lgamma(a) - (1.0 + a) * digamma(a)
end
function mgf(d::InverseGamma, t::Real)
a, b = d.shape, d.scale
a = d.shape
b = d.scale
(2 * (-b * t)^(a / 2)) / gamma(a) * besselk(a, sqrt(-4.0 * b * t))
end
function cf(d::InverseGamma, t::Real)
a, b = d.shape, d.scale
a = d.shape
b = d.scale
(2 * (-im * b * t)^(a / 2)) / gamma(a) * besselk(a, sqrt(-4.0 * im * b * t))
end
pdf(d::InverseGamma, x::Real) = exp(logpdf(d, x))
function logpdf(d::InverseGamma, x::Real)
a, b = d.shape, d.scale
-(a * log(b)) - lgamma(a) - ((a + 1.0) * log(x)) - 1.0/ (b * x)
a = d.shape
b = d.scale
a * log(b) - lgamma(a) - (a + 1.0) * log(x) - b / x
end
rand(d::InverseGamma) = 1.0 / rand(Gamma(d.shape, d.scale))
rand(d::InverseGamma) = 1.0 / rand(_inv(d))
function rand!(d::InverseGamma, A::Array{Float64})
rand!(Gamma(d.shape, d.scale), A)
rand!(_inv(d), A)
for i in 1:length(A)
A[i] = 1.0 / A[i]
end
......
......@@ -168,11 +168,11 @@ x = rand(Normal(2.0, 3.0), n)
p = posterior((2.0, pri), Normal, x)
@test isa(p, InverseGamma)
@test_approx_eq p.shape pri.shape + n / 2
@test_approx_eq rate(p) rate(pri) + sum(abs2(x - 2.0)) / 2
@test_approx_eq p.scale pri.scale + sum(abs2(x - 2.0)) / 2
p = posterior((2.0, pri), Normal, x, w)
@test isa(p, InverseGamma)
@test_approx_eq p.shape pri.shape + sum(w) / 2
@test_approx_eq rate(p) rate(pri) + dot(w, abs2(x - 2.0)) / 2
@test_approx_eq p.scale pri.scale + dot(w, abs2(x - 2.0)) / 2
......@@ -59,8 +59,8 @@ for d in [Arcsine(),
HyperGeometric(2.0, 2.0, 3.0),
InverseGaussian(1.0,1.0),
InverseGaussian(2.0,7.0),
InverseGamma(1.0,1.0),
InverseGamma(2.0,3.0),
InverseGamma(1.0, 1.0),
InverseGamma(2.0, 3.0),
# Kolmogorov(), # no quantile function
Laplace(0.0, 1.0),
Laplace(10.0, 1.0),
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment