Skip to content

Commit

Permalink
Check that PDSparseMat is defined in tests
Browse files Browse the repository at this point in the history
  • Loading branch information
ElOceanografo committed Nov 11, 2020
1 parent d349702 commit 26a53cb
Showing 1 changed file with 45 additions and 40 deletions.
85 changes: 45 additions & 40 deletions test/mvnormal.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,9 @@
# Tests on Multivariate Normal distributions

import PDMats: ScalMat, PDiagMat, PDMat, PDSparseMat
import PDMats: ScalMat, PDiagMat, PDMat
if isdefined(PDMats, :PDSparseMat)
import PDMats: PDSparseMat
end

using Distributions
using LinearAlgebra, Random, Test
Expand Down Expand Up @@ -178,45 +181,47 @@ end
end

##### Random sampling from MvNormalCanon with sparse precision matrix
@testset "Sparse MvNormalCanon random sampling" begin
n = 20
nsamp = 100_000
Random.seed!(1234)

J = sprandn(n, n, 0.25)
J = J'J + I
Σ = inv(Matrix(J))
J = PDSparseMat(J)
μ = zeros(n)

d_prec_sparse = MvNormalCanon(μ, J*μ, J)
d_prec_dense = MvNormalCanon(μ, J*μ, PDMat(Matrix(J)))
d_cov_dense = MvNormal(μ, PDMat(Symmetric(Σ)))

x_prec_sparse = rand(d_prec_sparse, nsamp)
x_prec_dense = rand(d_prec_dense, nsamp)
x_cov_dense = rand(d_cov_dense, nsamp)

dists = [d_prec_sparse, d_prec_dense, d_cov_dense]
samples = [x_prec_sparse, x_prec_dense, x_cov_dense]
tol = 1e-16
se = sqrt.(diag(Σ) ./ nsamp)
#=
The cholesky decomposition of sparse matrices is performed by `SuiteSparse.CHOLMOD`,
which returns a different decomposition than the `Base.LinearAlgebra` function (which uses
LAPACK). These different Cholesky routines produce different factorizations (since the
Cholesky factorization is not in general unique). As a result, the random samples from
an `MvNormalCanon` distribution with a sparse precision matrix are not in general
identical to those from an `MvNormalCanon` or `MvNormal`, even if the seeds are
identical. As a result, these tests only check for approximate statistical equality,
rather than strict numerical equality of the samples.
=#
for i in 1:3, j in 1:3
@test all(abs.(mean(samples[i]) .- μ) .< 2se)
loglik_ii = [logpdf(dists[i], samples[i][:, k]) for k in 1:100_000]
loglik_ji = [logpdf(dists[j], samples[i][:, k]) for k in 1:100_000]
# test average likelihood ratio between distribution i and sample j are small
@test mean((loglik_ii .- loglik_ji).^2) < tol
if isdefined(PDMats, :PDSparseMat)
@testset "Sparse MvNormalCanon random sampling" begin
n = 20
nsamp = 100_000
Random.seed!(1234)

J = sprandn(n, n, 0.25)
J = J'J + I
Σ = inv(Matrix(J))
J = PDSparseMat(J)
μ = zeros(n)

d_prec_sparse = MvNormalCanon(μ, J*μ, J)
d_prec_dense = MvNormalCanon(μ, J*μ, PDMat(Matrix(J)))
d_cov_dense = MvNormal(μ, PDMat(Symmetric(Σ)))

x_prec_sparse = rand(d_prec_sparse, nsamp)
x_prec_dense = rand(d_prec_dense, nsamp)
x_cov_dense = rand(d_cov_dense, nsamp)

dists = [d_prec_sparse, d_prec_dense, d_cov_dense]
samples = [x_prec_sparse, x_prec_dense, x_cov_dense]
tol = 1e-16
se = sqrt.(diag(Σ) ./ nsamp)
#=
The cholesky decomposition of sparse matrices is performed by `SuiteSparse.CHOLMOD`,
which returns a different decomposition than the `Base.LinearAlgebra` function (which uses
LAPACK). These different Cholesky routines produce different factorizations (since the
Cholesky factorization is not in general unique). As a result, the random samples from
an `MvNormalCanon` distribution with a sparse precision matrix are not in general
identical to those from an `MvNormalCanon` or `MvNormal`, even if the seeds are
identical. As a result, these tests only check for approximate statistical equality,
rather than strict numerical equality of the samples.
=#
for i in 1:3, j in 1:3
@test all(abs.(mean(samples[i]) .- μ) .< 2se)
loglik_ii = [logpdf(dists[i], samples[i][:, k]) for k in 1:100_000]
loglik_ji = [logpdf(dists[j], samples[i][:, k]) for k in 1:100_000]
# test average likelihood ratio between distribution i and sample j are small
@test mean((loglik_ii .- loglik_ji).^2) < tol
end
end
end

Expand Down

0 comments on commit 26a53cb

Please sign in to comment.