我有一个对称的 30 x 30 协方差矩阵(同时通过 isSymmetric 函数以及 Matrixcalc 库中的 is.symmetric.matrix 传递)。但是,我无法使用 rmvnorm() 从中采样,因为我收到以下错误:
chol.default(Sigma, ...) 中的错误: 14 阶的前导次要不是正数
重新创建矩阵 (V) 的代码如下。我看到其他人也有类似的问题,但他们中的大多数人意识到他们的矩阵只是勉强不对称。我的没有这个问题有修复吗?
V <- structure(c(29.7886188034286, 0.000135608020306483, 3.31533243601536e-05,
1.4149260534386e-08, 1.47707078131224e-10, 3.37580764203078e-07,
5.64221265795885e-13, 2.25017993462975e-15, 2.29195223482087e-16,
9.10148117311587e-22, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 0.000135608020306483, 29.7886188034286,
2.2426518115228e-08, 5.5775762335415e-12, 1.77927684962432e-14,
5.50469954608396e-09, 4.95218818529839e-17, 3.90846912833641e-11,
1.50295585574595e-17, 1.16206427894487e-20, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3.31533243601536e-05,
2.2426518115228e-08, 29.7886188034286, 0.00685475942462371, 2.11134810790327e-05,
0.0534401615284902, 5.72799389787849e-08, 4.18624110716291e-15,
1.11496083204703e-10, 2.41081638713732e-17, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.4149260534386e-08,
5.5775762335415e-12, 0.00685475942462371, 29.7886188034286, 0.0485441146570834,
0.00235235999633464, 0.000151744662859306, 1.92911292294946e-17,
1.2344146127685e-08, 9.27870135119539e-17, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.47707078131224e-10,
1.77927684962432e-14, 2.11134810790327e-05, 0.0485441146570834,
29.7886188034286, 4.12378769436744e-06, 0.077752722328188, 3.69961837029023e-20,
2.24120473537269e-09, 2.78180051438489e-18, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3.37580764203078e-07,
5.50469954608396e-09, 0.0534401615284902, 0.00235235999633464,
4.12378769436744e-06, 29.7886188034286, 1.53888416023551e-08,
2.27088584364582e-13, 1.90414928862951e-08, 1.27417002165266e-14,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5.64221265795885e-13,
4.95218818529839e-17, 5.72799389787849e-08, 0.000151744662859306,
0.077752722328188, 1.53888416023551e-08, 29.7886188034286, 1.87089710144973e-22,
6.07333453807195e-10, 2.7834529260172e-19, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.25017993462975e-15,
3.90846912833641e-11, 4.18624110716291e-15, 1.92911292294946e-17,
3.69961837029023e-20, 2.27088584364582e-13, 1.87089710144973e-22,
29.7886188034286, 1.09623746831197e-16, 7.73956371456588e-14,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.29195223482087e-16,
1.50295585574595e-17, 1.11496083204703e-10, 1.2344146127685e-08,
2.24120473537269e-09, 1.90414928862951e-08, 6.07333453807195e-10,
1.09623746831197e-16, 29.7886188034286, 9.14354495463089e-09,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 9.10148117311587e-22,
1.16206427894487e-20, 2.41081638713732e-17, 9.27870135119539e-17,
2.78180051438489e-18, 1.27417002165266e-14, 2.7834529260172e-19,
7.73956371456588e-14, 9.14354495463089e-09, 29.7886188034286,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1.1, 2.23435064143926e-11, 1.33547219558487e-12,
2.43247400092917e-19, 2.65083889750281e-23, 1.38463754600901e-16,
3.86793905062364e-28, 6.15198425841414e-33, 6.38251501624773e-35,
1.00647935339086e-45, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 2.23435064143926e-11, 1.1, 6.1108893902982e-19,
3.77982549438626e-26, 3.84651806542874e-31, 3.68169628548815e-20,
2.9797168689254e-36, 1.85606889272712e-24, 2.74456278291983e-37,
1.64074474950797e-43, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1.33547219558487e-12, 6.1108893902982e-19,
1.1, 5.70897210323724e-08, 5.41627098503697e-13, 3.46944376178704e-06,
3.98644533783302e-18, 2.1292624542656e-32, 1.51042770288842e-23,
7.06169554938132e-37, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 2.43247400092917e-19, 3.77982549438626e-26,
5.70897210323724e-08, 1.1, 2.86287682950135e-06, 6.72335136471411e-09,
2.79774037428088e-11, 4.52164040370437e-37, 1.85141095238044e-19,
1.04605642973501e-35, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 2.65083889750281e-23, 3.84651806542874e-31,
5.41627098503697e-13, 2.86287682950135e-06, 1.1, 2.06620454498461e-14,
7.3439529070571e-06, 1.66300890475797e-42, 6.10300580448813e-21,
9.40226930623985e-39, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1.38463754600901e-16, 3.68169628548815e-20,
3.46944376178704e-06, 6.72335136471411e-09, 2.06620454498461e-14,
1.1, 2.87734924851621e-19, 6.26572494547037e-29, 4.4053732445802e-19,
1.9725839084239e-31, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 3.86793905062364e-28, 2.9797168689254e-36,
3.98644533783302e-18, 2.79774037428088e-11, 7.3439529070571e-06,
2.87734924851621e-19, 1.1, 4.25285449747533e-47, 4.48162101889993e-22,
9.41344266929708e-41, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 6.15198425841414e-33, 1.85606889272712e-24,
2.1292624542656e-32, 4.52164040370437e-37, 1.66300890475797e-42,
6.26572494547037e-29, 4.25285449747533e-47, 1.1, 1.46012488822639e-35,
7.27802729314438e-30, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 6.38251501624773e-35, 2.74456278291983e-37,
1.51042770288842e-23, 1.85141095238044e-19, 6.10300580448813e-21,
4.4053732445802e-19, 4.48162101889993e-22, 1.46012488822639e-35,
1.1, 1.01580402447928e-19, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1.00647935339086e-45, 1.64074474950797e-43,
7.06169554938132e-37, 1.04605642973501e-35, 9.40226930623985e-39,
1.9725839084239e-31, 9.41344266929708e-41, 7.27802729314438e-30,
1.01580402447928e-19, 1.1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.1, 2.23435064143926e-11,
1.33547219558487e-12, 2.43247400092917e-19, 2.65083889750281e-23,
1.38463754600901e-16, 3.86793905062364e-28, 6.15198425841414e-33,
6.38251501624773e-35, 1.00647935339086e-45, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.23435064143926e-11,
1.1, 6.1108893902982e-19, 3.77982549438626e-26, 3.84651806542874e-31,
3.68169628548815e-20, 2.9797168689254e-36, 1.85606889272712e-24,
2.74456278291983e-37, 1.64074474950797e-43, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.33547219558487e-12,
6.1108893902982e-19, 1.1, 5.70897210323724e-08, 5.41627098503697e-13,
3.46944376178704e-06, 3.98644533783302e-18, 2.1292624542656e-32,
1.51042770288842e-23, 7.06169554938132e-37, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.43247400092917e-19,
3.77982549438626e-26, 5.70897210323724e-08, 1.1, 2.86287682950135e-06,
6.72335136471411e-09, 2.79774037428088e-11, 4.52164040370437e-37,
1.85141095238044e-19, 1.04605642973501e-35, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.65083889750281e-23,
3.84651806542874e-31, 5.41627098503697e-13, 2.86287682950135e-06,
1.1, 2.06620454498461e-14, 7.3439529070571e-06, 1.66300890475797e-42,
6.10300580448813e-21, 9.40226930623985e-39, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.38463754600901e-16,
3.68169628548815e-20, 3.46944376178704e-06, 6.72335136471411e-09,
2.06620454498461e-14, 1.1, 2.87734924851621e-19, 6.26572494547037e-29,
4.4053732445802e-19, 1.9725839084239e-31, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3.86793905062364e-28,
2.9797168689254e-36, 3.98644533783302e-18, 2.79774037428088e-11,
7.3439529070571e-06, 2.87734924851621e-19, 1.1, 4.25285449747533e-47,
4.48162101889993e-22, 9.41344266929708e-41, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6.15198425841414e-33,
1.85606889272712e-24, 2.1292624542656e-32, 4.52164040370437e-37,
1.66300890475797e-42, 6.26572494547037e-29, 4.25285449747533e-47,
1.1, 1.46012488822639e-35, 7.27802729314438e-30, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6.38251501624773e-35,
2.74456278291983e-37, 1.51042770288842e-23, 1.85141095238044e-19,
6.10300580448813e-21, 4.4053732445802e-19, 4.48162101889993e-22,
1.46012488822639e-35, 1.1, 1.01580402447928e-19, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.00647935339086e-45,
1.64074474950797e-43, 7.06169554938132e-37, 1.04605642973501e-35,
9.40226930623985e-39, 1.9725839084239e-31, 9.41344266929708e-41,
7.27802729314438e-30, 1.01580402447928e-19, 1.1), dim = c(30L,
30L))
rmvnorm(1, sigma = V)
你的问题不在于缺乏对称性,而在于缺乏正定性。这也不是矩阵“几乎”正定的常见情况(即最小特征值接近零,无论是正值还是负值);最小的特征值显然是负的(-8.9,相比之下主要特征值是 37.4)。
eigen(V)$values
[1] 37.424265 29.861301 29.823050 29.788700 29.788619 29.788619 29.788483
[8] 29.788059 29.735071 29.697642 3.502381 1.100006 1.100006 1.100002
[15] 1.100002 1.100000 1.100000 1.100000 1.100000 1.100000 1.100000
[22] 1.100000 1.100000 1.100000 1.100000 1.099997 1.099997 1.099992
[29] 1.099992 -8.899997
mvtnorm::rmvnorm
将让您使用此协方差矩阵对 MVN 值进行采样,但在警告您之后(“西格玛在数值上不是正半定的”)
它将与负特征值对应的所有分量截断为零。由您决定这对您的应用程序是否有意义(!!) 无论如何,R 生态系统中有很多
rmvnorm()
函数 — 这些只是我在系统上安装的包含
rmvnorm()
函数的软件包。他们都会以不同的方式处理非正定协方差矩阵......unique(help.search("rmvnorm", agrep = FALSE)$matches$Package)
## [1] "CDM" "dae" "dqrng" "geoR"
## [5] "ks" "mhsmm" "miceadds" "mixtools"
## [9] "mvtnorm" "rethinking" "Rfast" "spam"
## [13] "splus2R" "stockassessment"