From 26a53cb7c776028bb79797518aa3676b7704313b Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Wed, 11 Nov 2020 09:19:02 -0800 Subject: [PATCH] Check that PDSparseMat is defined in tests --- test/mvnormal.jl | 85 +++++++++++++++++++++++++----------------------- 1 file changed, 45 insertions(+), 40 deletions(-) diff --git a/test/mvnormal.jl b/test/mvnormal.jl index 9e34ec181..be3d590e0 100644 --- a/test/mvnormal.jl +++ b/test/mvnormal.jl @@ -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 @@ -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