Ce serveur Gitlab sera éteint le 30 juin 2020, pensez à migrer vos projets vers les serveurs gitlab-research.centralesupelec.fr et gitlab-student.centralesupelec.fr !

Commit d5da72d0 authored by Dahua Lin's avatar Dahua Lin

make canonical form of multivariate normal into a full-fledged distribution

parent cd475a2c
......@@ -48,6 +48,8 @@ export
Exponential,
FDist,
Gamma,
GenericMvNormal,
GenericMvNormalCanon,
Geometric,
Gumbel,
HyperGeometric,
......@@ -265,7 +267,8 @@ include(joinpath("univariate", "weibull.jl"))
# Multivariate distributions
include(joinpath("multivariate", "dirichlet.jl"))
include(joinpath("multivariate", "multinomial.jl"))
include(joinpath("multivariate", "multivariatenormal.jl"))
include(joinpath("multivariate", "mvnormal.jl"))
include(joinpath("multivariate", "mvnormalcanon.jl"))
include(joinpath("multivariate", "vonmisesfisher.jl"))
# Matrix distributions
......
......@@ -2,7 +2,7 @@
#### Generic MvNormal -- Generic MvNormal (Σ is known)
function posterior_canon(prior::GenericMultivariateNormal, ss::GenericMvNormalKnownSigmaStats)
function posterior_canon(prior::GenericMvNormal, ss::GenericMvNormalKnownSigmaStats)
invΣ0 = inv(prior.Σ)
μ0 = prior.μ
invΣp = add_scal(invΣ0, ss.invΣ, ss.tw)
......@@ -10,46 +10,46 @@ function posterior_canon(prior::GenericMultivariateNormal, ss::GenericMvNormalKn
return GenericMvNormalCanon(h, invΣp)
end
function posterior_canon{Pri<:GenericMultivariateNormal,Cov<:AbstractPDMat}(
function posterior_canon{Pri<:GenericMvNormal,Cov<:AbstractPDMat}(
prior::(Pri, Cov),
G::Type{GenericMultivariateNormal{Cov}},
G::Type{GenericMvNormal{Cov}},
x::Matrix)
μpri::Pri, Σ::Cov = prior
posterior_canon(μpri, suffstats(GenericMvNormalKnownSigma{Cov}(Σ), x))
end
function posterior_canon{Pri<:GenericMultivariateNormal,Cov<:AbstractPDMat}(
function posterior_canon{Pri<:GenericMvNormal,Cov<:AbstractPDMat}(
prior::(Pri, Cov),
G::Type{GenericMultivariateNormal{Cov}},
G::Type{GenericMvNormal{Cov}},
x::Matrix, w::Array{Float64})
μpri::Pri, Σ::Cov = prior
posterior_canon(μpri, suffstats(GenericMvNormalKnownSigma{Cov}(Σ), x, w))
end
function posterior{Pri<:GenericMultivariateNormal,Cov<:AbstractPDMat}(
function posterior{Pri<:GenericMvNormal,Cov<:AbstractPDMat}(
prior::(Pri, Cov),
G::Type{GenericMultivariateNormal{Cov}},
G::Type{GenericMvNormal{Cov}},
x::Matrix)
convert(Pri, posterior_canon(prior, G, x))
end
function posterior{Pri<:GenericMultivariateNormal,Cov<:AbstractPDMat}(
function posterior{Pri<:GenericMvNormal,Cov<:AbstractPDMat}(
prior::(Pri, Cov),
G::Type{GenericMultivariateNormal{Cov}},
G::Type{GenericMvNormal{Cov}},
x::Matrix, w::Array{Float64})
convert(Pri, posterior_canon(prior, G, x, w))
end
function complete{Pri<:GenericMultivariateNormal,Cov<:AbstractPDMat}(
G::Type{GenericMultivariateNormal{Cov}},
function complete{Pri<:GenericMvNormal,Cov<:AbstractPDMat}(
G::Type{GenericMvNormal{Cov}},
pri::(Pri, Cov),
μ::Vector{Float64})
GenericMultivariateNormal(μ, pri[2]::Cov)
GenericMvNormal(μ, pri[2]::Cov)
end
......@@ -64,7 +64,7 @@ function posterior(pri::(MvNormal, Matrix{Float64}), G::Type{MvNormal}, args...)
end
function complete(G::Type{MvNormal}, pri::(MvNormal, Matrix{Float64}), μ::Vector{Float64})
GenericMultivariateNormal(μ, PDMat(pri[2]))
GenericMvNormal(μ, PDMat(pri[2]))
end
......
# Canonical form of multivariate normal
## generic types
immutable GenericMvNormalCanon{Prec<:AbstractPDMat} <: AbstractMvNormal
dim::Int # dimension of sample space
zeromean::Bool # whether the mean vector is zero
μ::Vector{Float64} # the mean vector
h::Vector{Float64} # potential vector, i.e. inv(Σ) * μ
J::Prec # precision matrix, i.e. inv(Σ)
function GenericMvNormalCanon(μ::Vector{Float64}, h::Vector{Float64}, J::Prec, zmean::Bool)
d = dim(J)
length(μ) == length(h) == d || throw(ArgumentError("Inconsistent argument dimensions."))
new(d, zmean, μ, h, J)
end
function GenericMvNormalCanon(J::Prec)
d = dim(J)
new(d, true, zeros(d), zeros(d), J)
end
end
function GenericMvNormalCanon{P<:AbstractPDMat}(μ::Vector{Float64}, h::Vector{Float64}, J::P)
GenericMvNormalCanon{P}(μ, h, J, allzeros(μ))
end
function GenericMvNormalCanon{P<:AbstractPDMat}(h::Vector{Float64}, J::P, zmean::Bool)
μ = zmean ? zeros(length(h)) : J \ h
GenericMvNormalCanon{P}(μ, h, J, zmean)
end
function GenericMvNormalCanon{P<:AbstractPDMat}(h::Vector{Float64}, J::P)
GenericMvNormalCanon(h, J, allzeros(h))
end
GenericMvNormalCanon{P<:AbstractPDMat}(J::P) = GenericMvNormalCanon{P}(J)
## type aliases and convenient constructors
typealias MvNormalCanon GenericMvNormalCanon{PDMat}
typealias DiagNormalCanon GenericMvNormalCanon{PDiagMat}
typealias IsoNormalCanon GenericMvNormalCanon{ScalMat}
MvNormalCanon(J::PDMat) = GenericMvNormalCanon(J)
MvNormalCanon(J::Matrix{Float64}) = GenericMvNormalCanon(PDMat(J))
MvNormalCanon(h::Vector{Float64}, J::PDMat) = GenericMvNormalCanon(h, J)
MvNormalCanon(h::Vector{Float64}, J::Matrix{Float64}) = GenericMvNormalCanon(h, PDMat(J))
DiagNormalCanon(J::PDiagMat) = GenericMvNormalCanon(J)
DiagNormalCanon(J::Vector{Float64}) = GenericMvNormalCanon(PDiagMat(J))
DiagNormalCanon(h::Vector{Float64}, J::PDiagMat) = GenericMvNormalCanon(h, J)
DiagNormalCanon(h::Vector{Float64}, J::Vector{Float64}) = GenericMvNormalCanon(h, PDiagMat(J))
IsoNormalCanon(J::ScalMat) = GenericMvNormalCanon(J)
IsoNormalCanon(d::Integer, prec::Real) = GenericMvNormalCanon(ScalMat(int(d), float64(prec)))
IsoNormalCanon(h::Vector{Float64}, J::ScalMat) = GenericMvNormalCanon(h, J)
IsoNormalCanon(h::Vector{Float64}, prec::Real) = GenericMvNormalCanon(h, ScalMat(length(h), float64(prec)))
# conversion between conventional form and canonical form
function Base.convert{C<:AbstractPDMat}(D::Type{GenericMvNormal{C}}, cf::GenericMvNormalCanon{C})
GenericMvNormal{C}(cf.dim, cf.zeromean, cf.μ, inv(cf.J))
end
function Base.convert{C<:AbstractPDMat}(D::Type{GenericMvNormalCanon{C}}, d::GenericMvNormal{C})
J::C = inv(d.Σ)
h::Vector{Float64} = J * d.μ
GenericMvNormalCanon{C}(d.μ, h, J, d.zeromean)
end
canonform{C<:AbstractPDMat}(d::GenericMvNormal{C}) = convert(GenericMvNormalCanon{C}, d)
# Basic statistics
dim(d::GenericMvNormalCanon) = d.dim
mean(d::GenericMvNormalCanon) = d.μ
mode(d::GenericMvNormalCanon) = d.μ
modes(d::GenericMvNormalCanon) = [mode(d)]
var(d::GenericMvNormalCanon) = diag(inv(d.J))
cov(d::GenericMvNormalCanon) = full(inv(d.J))
invcov(d::GenericMvNormalCanon) = full(d.J)
logdet_cov(d::GenericMvNormalCanon) = -logdet(d.J)
entropy(d::GenericMvNormalCanon) = 0.5 * (dim(d) * (float64(log2π) + 1.0) - logdet(d.J))
# PDF evaluation
function sqmahal(d::GenericMvNormalCanon, x::Vector{Float64})
z::Vector{Float64} = d.zeromean ? x : x - d.μ
quad(d.J, z)
end
function sqmahal!(r::Array{Float64}, d::GenericMvNormalCanon, x::Matrix{Float64})
if !(size(x, 1) == dim(d) && size(x, 2) == length(r))
throw(ArgumentError("Inconsistent argument dimensions."))
end
z::Matrix{Float64} = d.zeromean ? x : bsubtract(x, d.μ, 1)
quad!(r, d.J, z)
end
# Sampling (for GenericMvNormal)
# helper routines (unwhiten w.r.t. inv(J))
unwhiten_winv!(J::ScalMat, z::VecOrMat{Float64}) = whiten!(J, z)
unwhiten_winv!(J::PDiagMat, z::VecOrMat{Float64}) = whiten!(J, z)
unwhiten_winv!(J::PDMat, x::VecOrMat{Float64}) = (Base.LinAlg.LAPACK.trtrs!('U', 'N', 'N', J.chol.UL, x); x)
# sampling functions
function rand!(d::GenericMvNormalCanon, x::Vector{Float64})
unwhiten_winv!(d.J, randn!(x))
if !d.zeromean
add!(x, d.μ)
end
x
end
function rand!(d::GenericMvNormalCanon, x::Matrix{Float64})
unwhiten_winv!(d.J, randn!(x))
if !d.zeromean
badd!(x, d.μ, 1)
end
x
end
......@@ -11,6 +11,15 @@ function allfinite{T<:Real}(x::Array{T})
return true
end
function allzeros{T<:Real}(x::Array{T})
for i = 1 : length(x)
if !(x == zero(T))
return false
end
end
return true
end
function isprobvec(p::Vector{Float64})
s = 0.
for i = 1:length(p)
......
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