Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixed posterior for NormalWishart. #28

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Fixed posterior for NormalWishart. #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
13ba68b
Fixed posterior for NormalWishart.
kaselby Apr 17, 2018
a323067
Updated to change parametrization.
kaselby May 2, 2018
c14f301
updated pdf function
kaselby May 5, 2018
57f52b1
whoops, thought I included these changes in the earlier commit.
kaselby May 5, 2018
098841d
Fixed.
kaselby May 5, 2018
5e5e4d6
updated to return to original parameterization.
kaselby Jun 14, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions src/mvnormal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,9 +106,8 @@ function posterior_canon(prior::NormalWishart, ss::MvNormalStats)
nu = nu0 + ss.tw
mu = (kappa0.*mu0 + ss.s) ./ kappa

Lam0 = TC0[:U]'*TC0[:U]
z = prior.zeromean ? ss.m : ss.m - mu0
Lam = Lam0 + ss.s2 + kappa0*ss.tw/kappa*(z*z')
Lam = Symmetric(inv(inv(TC0) + ss.s2 + kappa0*ss.tw/kappa*(z*z')))

return NormalWishart(mu, kappa, cholfact(Lam), nu)
end
Expand Down
20 changes: 11 additions & 9 deletions src/normalwishart.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ struct NormalWishart{T<:Real} <: ContinuousMultivariateDistribution
zeromean::Bool
mu::Vector{T}
kappa::T
Tchol::Cholesky{T} # Precision matrix (well, sqrt of one)
Tchol::Cholesky{T} #Prior inverse scale matrix
nu::T

function NormalWishart{T}(mu::Vector{T}, kappa::T,
Expand Down Expand Up @@ -40,7 +40,7 @@ end
function insupport(::Type{NormalWishart}, x::Vector{T}, Lam::Matrix{T}) where T<:Real
return (all(isfinite(x)) &&
size(Lam, 1) == size(Lam, 2) &&
isApproxSymmmetric(Lam) &&
#isApproxSymmmetric(Lam) &&
size(Lam, 1) == length(x) &&
hasCholesky(Lam))
end
Expand All @@ -52,7 +52,7 @@ function logpdf(nw::NormalWishart, x::Vector{T}, Lam::Matrix{T}) where T<:Real
if !insupport(NormalWishart, x, Lam)
return -Inf
else
p = length(x)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

did you mean to remove this? seems like it's needed below...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

alternatively, you could use nw.dim

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this was an accident

p=length(x)

nu = nw.nu
kappa = nw.kappa
Expand All @@ -62,10 +62,10 @@ function logpdf(nw::NormalWishart, x::Vector{T}, Lam::Matrix{T}) where T<:Real
hp = 0.5 * p

# Normalization
logp = hp*(log(kappa) - Float64(log2π))
logp = hp*(log(kappa) - Float64(log(2*π)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why get rid of log2pi? It's a StatsFuns constant.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My bad, I think it was giving me errors when I tried to run it. Maybe I just wasn't importing the required package or something.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I think it's necessary to add StatsFuns (which also provides logmvgamma which is needed elsewhere).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

logmvgamma is from Distributions. maybe the old lpgamma is a version from StatsFuns?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think Distributions gets it from StatsFuns: https://github.com/JuliaStats/StatsFuns.jl/blob/e66dd973650c375bc1739c820e5b96bb5bd000a8/src/misc.jl#L3-L15

logp -= hnu * logdet(Tchol)
logp -= hnu * p * log(2.)
logp -= lpgamma(p, hnu)
logp -= logmvgamma(p, hnu)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm a little disturbed that this hadn't been fixed until now......

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There were other areas too where the current logpdf function crashes. I think somewhere in insupport there is a call to a function that doesn't exist. I didn't fix that in this pull request because it was a separate issue, and I forgot to make a new issue for it.


# Wishart (MvNormal contributes 0.5 as well)
logp += (hnu - hp) * logdet(Lam)
Expand All @@ -76,13 +76,15 @@ function logpdf(nw::NormalWishart, x::Vector{T}, Lam::Matrix{T}) where T<:Real
logp -= 0.5 * kappa * dot(z, Lam * z)

return logp

end
end

function rand(nw::NormalWishart)
Lam = rand(Wishart(nw.nu, nw.Tchol))
Lsym = PDMat(Symmetric(inv(Lam) ./ nw.kappa))
mu = rand(MvNormal(nw.mu, Lsym))
# forcibly symmetrize to pass checks in Distributions.Wishart
V = PDMat(Symmetric(inv(nw.Tchol)))
Lam = rand(Wishart(nw.nu, V))
# forcibly symmetrize to pass checks in Distributions.MvNormal
Σsym = PDMat(Symmetric(inv(Lam) ./ nw.kappa))
mu = rand(MvNormal(nw.mu, Σsym))
return (mu, Lam)
end