MvNormal Ошибка с симметричной и неотрицательно определённой матрицей

Question

сводки моей проблемы в том, что я пытаюсь повторить функцию Matlab:

mvnrnd(mu', sigma, 200) 

в Джулию с помощью:

rand(MvNormal(mu, sigma), 200)' 

и результата представляет собой матрицу 200 x 7, по существу генерирующую 200 случайных временных рядов данных.

Matlab works, Julia нет.

Мои входные матрицы:

mu = [0.15; 0.03; 0.06; 0.04; 0.1; 0.02; 0.12] 

sigma = [0.0035 -0.0038 0.0020 0.0017 -0.0006 -0.0028 0.0009; 
    -0.0038 0.0046 -0.0011 0.0001 0.0003 0.0054 -0.0024; 
    0.0020 -0.0011 0.0041 0.0068 -0.0004 0.0047 -0.0036; 
    0.0017 0.0001 0.0068 0.0125 0.0002 0.0109 -0.0078; 
    -0.0006 0.0003 -0.0004 0.0002 0.0025 -0.0004 -0.0007; 
    -0.0028 0.0054 0.0047 0.0109 -0.0004 0.0159 -0.0093; 
    0.0009 -0.0024 -0.0036 -0.0078 -0.0007 -0.0093 0.0061] 

Использование Distributions.jl, запустив линию:

MvNormal(sigma) 

производит ошибку:

ERROR: LoadError: Base.LinAlg.PosDefException(4) 

Матрица сигма является симметричным, но только положительный полуопределенный:

issym(sigma) #symmetrical 
> true 
isposdef(sigma) #positive definite 
> false 

using LinearOperators 
check_positive_definite(sigma) #check for positive (semi-)definite 
> true 

Matlab дает те же результаты для этих тестов, однако Matlab способен генерировать матрицу случайного возврата 200x7.

Может кто-нибудь посоветует, что я могу сделать, чтобы заставить его работать в Джулии? Или где проблема?

Спасибо.

Answer 1

Вопрос заключается в том, что ковариационная матрица является неопределенной. См.

julia> eigvals(sigma) 
7-element Array{Float64,1}: 
-3.52259e-5 
-2.42008e-5 
    2.35508e-7 
    7.08269e-5 
    0.00290538 
    0.0118957 
    0.0343873 

так что это не ковариационная матрица. Возможно, это произошло из-за округления, поэтому, если у вас есть доступ к некруглым данным, вы можете попробовать это вместо этого. Я просто попытался, и у меня также появилась ошибка в Matlab. Однако, в отличие от Джулии, Matlab позволяет матрице быть положительной полу определенно.

Способ сделать эту работу - добавить диагональную матрицу к исходной матрице, а затем ввести ее в MvNormal. То есть

julia> MvNormal(randn(7), sigma - minimum(eigvals(Symmetric(sigma)))*I) 
Distributions.MvNormal{PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}}(
dim: 7 
μ: [0.889004,-0.768551,1.78569,0.130445,0.589029,0.529418,-0.258474] 
Σ: 7x7 Array{Float64,2}: 
    0.00353523 -0.0038  0.002  0.0017  -0.0006  -0.0028  0.0009  
-0.0038  0.00463523 -0.0011  0.0001  0.0003  0.0054  -0.0024  
    0.002  -0.0011  0.00413523 0.0068  -0.0004  0.0047  -0.0036  
    0.0017  0.0001  0.0068  0.0125352 0.0002  0.0109  -0.0078  
-0.0006  0.0003  -0.0004  0.0002  0.00253523 -0.0004  -0.0007  
-0.0028  0.0054  0.0047  0.0109  -0.0004  0.0159352 -0.0093  
    0.0009  -0.0024  -0.0036  -0.0078  -0.0007  -0.0093  0.00613523 
) 

«Ковариационная» матрица, конечно, не то же самое, но она очень близка.