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 08aa0e9b authored by Dahua Lin's avatar Dahua Lin

tests for mvnormalcanon

parent d5da72d0
......@@ -37,6 +37,7 @@ export
Chisq,
Cosine,
DiagNormal,
DiagNormalCanon,
Dirichlet,
DiscreteUniform,
DoubleExponential,
......@@ -56,7 +57,8 @@ export
InverseWishart,
InverseGamma,
InverseGaussian,
IsoNormal,
IsoNormal,
IsoNormalCanon,
Kolmogorov,
KSDist,
KSOneSided,
......@@ -68,6 +70,7 @@ export
Multinomial,
MultivariateNormal,
MvNormal,
MvNormalCanon,
MvNormalKnownSigma,
NegativeBinomial,
NoncentralBeta,
......@@ -112,6 +115,7 @@ export
freecumulant, # free cumulants of distribution
gmvnormal, # a generic function to construct multivariate normal distributions
insupport, # predicate, is x in the support of the distribution?
invcov, # get the inversed covariance
invlogccdf, # complementary quantile based on log probability
invlogcdf, # quantile based on log probability
isplatykurtic, # Is excess kurtosis > 0.0?
......
......@@ -26,7 +26,7 @@ function GenericMvNormalCanon{P<:AbstractPDMat}(μ::Vector{Float64}, h::Vector{F
end
function GenericMvNormalCanon{P<:AbstractPDMat}(h::Vector{Float64}, J::P, zmean::Bool)
μ = zmean ? zeros(length(h)) : J \ h
μ = zmean ? zeros(length(h)) : (J \ h)
GenericMvNormalCanon{P}(μ, h, J, zmean)
end
......@@ -34,7 +34,10 @@ function GenericMvNormalCanon{P<:AbstractPDMat}(h::Vector{Float64}, J::P)
GenericMvNormalCanon(h, J, allzeros(h))
end
GenericMvNormalCanon{P<:AbstractPDMat}(J::P) = GenericMvNormalCanon{P}(J)
function GenericMvNormalCanon{P<:AbstractPDMat}(J::P)
d = dim(J)
GenericMvNormalCanon{P}(zeros(d), zeros(d), J, true)
end
## type aliases and convenient constructors
......@@ -106,14 +109,6 @@ 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
......
......@@ -5,13 +5,19 @@ const n_samples = 5_000_001
mu = [1.0, 2.0, 3.0]
C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]
h = mu
J = C
for d in [
Dirichlet(3, 2.0),
Dirichlet([2.0, 1.0, 3.0]),
MultivariateNormal(mu, 2.0),
MultivariateNormal(mu, [1.5, 2.0, 2.5]),
MultivariateNormal(mu, C)]
IsoNormal(mu, 2.0),
DiagNormal(mu, [1.5, 2.0, 2.5]),
MvNormal(mu, C),
IsoNormalCanon(h, 2.0),
DiagNormalCanon(h, [1.5, 2.0, 1.2]),
MvNormalCanon(h, J)]
println(d)
dmean = mean(d)
......
# Canonical form of multivariate normal
import NumericExtensions
import NumericExtensions.ScalMat
import NumericExtensions.PDiagMat
import NumericExtensions.PDMat
using Distributions
using Base.Test
##### construction, basic properties, and evaluation
h = [1., 2., 3.]
dv = [1.2, 3.4, 2.6]
J = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]
x1 = [3.2, 1.8, 2.4]
x = rand(3, 100)
# SGauss
gs = IsoNormalCanon(h, 2.0)
@test isa(gs, IsoNormalCanon)
@test dim(gs) == 3
@test mean(gs) == mode(gs) == h / 2.0
@test !gs.zeromean
@test invcov(gs) == diagm(fill(2.0, 3))
@test cov(gs) == diagm(fill(0.5, 3))
@test var(gs) == diag(cov(gs))
@test_approx_eq entropy(gs) 0.5 * logdet(2π * e * cov(gs))
gsz = IsoNormalCanon(3, 2.0)
@test isa(gsz, IsoNormalCanon)
@test dim(gsz) == 3
@test mean(gsz) == zeros(3)
@test gsz.zeromean
# DGauss
gd = DiagNormalCanon(h, dv)
@test isa(gd, DiagNormalCanon)
@test dim(gd) == 3
@test_approx_eq mean(gd) h ./ dv
@test !gd.zeromean
@test invcov(gd) == diagm(dv)
@test_approx_eq cov(gd) diagm(1.0 / dv)
@test_approx_eq var(gd) diag(cov(gd))
@test_approx_eq entropy(gd) 0.5 * logdet(2π * e * cov(gd))
gdz = DiagNormalCanon(dv)
@test isa(gdz, DiagNormalCanon)
@test dim(gdz) == 3
@test mean(gdz) == zeros(3)
@test gdz.zeromean
# Gauss
gf = MvNormalCanon(h, J)
@test isa(gf, MvNormalCanon)
@test dim(gf) == 3
@test_approx_eq mean(gf) J \ h
@test !gf.zeromean
@test invcov(gf) == J
@test_approx_eq cov(gf) inv(J)
@test_approx_eq var(gf) diag(cov(gf))
@test_approx_eq entropy(gf) 0.5 * logdet(2π * e * cov(gf))
gfz = MvNormalCanon(J)
@test isa(gfz, MvNormalCanon)
@test dim(gfz) == 3
@test mean(gfz) == zeros(3)
@test gfz.zeromean
# conversion
us = convert(IsoNormal, gs)
gs2 = convert(IsoNormalCanon, us)
@test isa(us, IsoNormal)
@test isa(gs2, IsoNormalCanon)
@test dim(gs2) == dim(gs)
@test_approx_eq gs2.h gs.h
@test_approx_eq gs2.J.value gs.J.value
ud = convert(DiagNormal, gd)
gd2 = convert(DiagNormalCanon, ud)
@test isa(ud, DiagNormal)
@test isa(gd2, DiagNormalCanon)
@test dim(gd2) == dim(gd)
@test_approx_eq gd2.h gd.h
@test_approx_eq gd2.J.diag gd.J.diag
uf = convert(MvNormal, gf)
gf2 = convert(MvNormalCanon, uf)
@test isa(uf, MvNormal)
@test isa(gf2, MvNormalCanon)
@test dim(gf2) == dim(gf)
@test_approx_eq gf2.h gf.h
@test_approx_eq gf2.J.mat gf.J.mat
# logpdf evaluation
@test_approx_eq logpdf(gs, x1) logpdf(us, x1)
@test_approx_eq logpdf(gs, x) logpdf(us, x)
@test_approx_eq logpdf(gsz, x1) logpdf(convert(IsoNormal, gsz), x1)
@test_approx_eq logpdf(gsz, x) logpdf(convert(IsoNormal, gsz), x)
@test_approx_eq logpdf(gd, x1) logpdf(ud, x1)
@test_approx_eq logpdf(gd, x) logpdf(ud, x)
@test_approx_eq logpdf(gdz, x1) logpdf(convert(DiagNormal, gdz), x1)
@test_approx_eq logpdf(gdz, x) logpdf(convert(DiagNormal, gdz), x)
@test_approx_eq logpdf(gf, x1) logpdf(uf, x1)
@test_approx_eq logpdf(gf, x) logpdf(uf, x)
@test_approx_eq logpdf(gfz, x1) logpdf(convert(MvNormal, gfz), x1)
@test_approx_eq logpdf(gfz, x) logpdf(convert(MvNormal, gfz), x)
# sampling
x = rand(gs)
@test isa(x, Vector{Float64})
@test length(x) == dim(gs)
x = rand(gd)
@test isa(x, Vector{Float64})
@test length(x) == dim(gd)
x = rand(gf)
@test isa(x, Vector{Float64})
@test length(x) == dim(gf)
x = rand(gsz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gsz)
x = rand(gdz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gdz)
x = rand(gfz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gfz)
n = 10
x = rand(gs, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gs), n)
x = rand(gd, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gd), n)
x = rand(gf, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gf), n)
x = rand(gsz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gsz), n)
x = rand(gdz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gdz), n)
x = rand(gfz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gfz), n)
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