我的问题的总结是我正在尝试复制 Matlab 函数:
mvnrnd(mu', sigma, 200)
进入朱莉娅使用:
rand( MvNormal(mu, sigma), 200)'
结果是一个 200 x 7 的矩阵,基本上生成 200 个随机返回时间序列数据。
Matlab有效,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 的随机返回样本矩阵。
有人可以建议我能做些什么来让它在朱莉娅工作吗?或者问题出在哪里?
谢谢。
问题是协方差矩阵是不确定的。看
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 中也遇到了错误。然而,与Julia相比,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
)
"协方差"矩阵当然不再相同,但它非常接近。