How do I fit an MvNormal to a matrix with Distributions.jl?

General Usage
, ,
6

You have more variables than observations so you’re hitting a problem with Distributions.jl: it only supports covariance matrices that are positive definite.

In your case the observation vectors are obviously linearly dependent so the covariance matrix is positive semidefinite. This is easy to reproduce:

julia> MvNormal(cov(randn(3, 5)))
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.

Several packages have been created to work around that limitation. You can try GitHub - invenia/PDMatsExtras.jl: Extra Positive (Semi-)Definite Matricies which allows you to wrap your matrix in a PSDMat type:

julia> MvNormal(PSDMat(cov(randn(3, 5))))
MvNormal{Float64, PSDMat{Float64, Matrix{Float64}}, FillArrays.Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}}(
dim: 5
μ: 5-element Zeros{Float64}
Σ: [1.8178793715225665 -0.10575355536032809 … 1.7031492738916045 -0.44842050805508626; -0.10575355536032809 0.055479378469435375 … -0.12077962629136424 0.10681936920348316; … ; 1.7031492738916045 -0.12077962629136424 … 1.6052066188918248 -0.4556363541776913; -0.44842050805508626 0.10681936920348316 … -0.4556363541776913 0.2427467729480703]
)

although as I understand that’s quite a hack: PSDMat is a subtype of AbstractPDMat so it pretends the matrix is positive definite. I don’t know if that can cause problems elsewhere…

Depending on your use case another option is to use the Gaussian type from GitHub - mschauer/GaussianDistributions.jl: Gaussian distributions as state variables and for uncertainty quantification with Unitful

2 Likes