2

我一直在尝试在 R 中重现类似 Cholesky 的协方差分解 - 就像在 Matlab 中使用cholcov(). 示例取自https://uk.mathworks.com/help/stats/cholcov.html。

作为示例的原始cholcov()函数的结果：

T =
   -0.2113    0.7887   -0.5774         0
    0.7887   -0.2113   -0.5774         0
    1.1547    1.1547    1.1547    1.7321

我试图T在 R 中复制它。我试过：

C1 <- cbind(c(2,1,1,2), c(1,2,1,2), c(1,1,2,2), c(2,2,2,3))
T1 <- chol(C1)
C2 <- t(T1) %*% T1

我的结果：

         [,1]      [,2]      [,3]         [,4]
[1,] 1.414214 0.7071068 0.7071068 1.414214e+00
[2,] 0.000000 1.2247449 0.4082483 8.164966e-01
[3,] 0.000000 0.0000000 1.1547005 5.773503e-01
[4,] 0.000000 0.0000000 0.0000000 1.290478e-08

C2恢复C1，但T1与 MATLAB 的解决方案完全不同。然后我想也许这将是协方差矩阵的 Cholesky 组合：

T1 <- chol(cov(C1))

但我明白了

          [,1]      [,2]      [,3]         [,4]
[1,] 0.5773503 0.0000000 0.0000000 2.886751e-01
[2,] 0.0000000 0.5773503 0.0000000 2.886751e-01
[3,] 0.0000000 0.0000000 0.5773503 2.886751e-01
[4,] 0.0000000 0.0000000 0.0000000 3.725290e-09

这也不对。

谁能给我一个提示如何cholcov()在 Matlab 中计算以便我可以在 R 中复制它？

4

1 回答 1

3

chol在这种情况下，您实际上是在滥用 R 功能。cholcovMATLAB 中的函数是复合函数。

  • 如果协方差为正，则进行 Cholesky 分解，返回全秩上三角 Cholesky 因子；
  • 如果协方差是半正定的，它会进行特征分解，返回一个矩形矩阵。

另一方面，cholfrom R 只进行 Choleksy 分解。您给出的示例C1属于第二种情况。因此，我们应该求助于eigenR 中的函数。

E <- eigen(C1, symmetric = TRUE)
#$values
#[1] 7.000000e+00 1.000000e+00 1.000000e+00 2.975357e-17
#
#$vectors
#           [,1]          [,2]          [,3]       [,4]
#[1,] -0.4364358  0.000000e+00  8.164966e-01 -0.3779645
#[2,] -0.4364358 -7.071068e-01 -4.082483e-01 -0.3779645
#[3,] -0.4364358  7.071068e-01 -4.082483e-01 -0.3779645
#[4,] -0.6546537  8.967707e-16 -2.410452e-16  0.7559289

V <- E$vectors
D <- sqrt(E$values)  ## root eigen values

由于数值秩为 3，我们删除最后一个特征值和特征向量：

V1 <- V[, 1:3]
D1 <- D[1:3]

在此处输入图像描述

因此，您想要的因素是：

R <- D1 * t(V1)  ## diag(D1) %*% t(V1)
#           [,1]       [,2]       [,3]          [,4]
#[1,] -1.1547005 -1.1547005 -1.1547005 -1.732051e+00
#[2,]  0.0000000 -0.7071068  0.7071068  8.967707e-16
#[3,]  0.8164966 -0.4082483 -0.4082483 -2.410452e-16

我们可以验证：

crossprod(R)  ## t(R) %*% R

#     [,1] [,2] [,3] [,4]
#[1,]    2    1    1    2
#[2,]    1    2    1    2
#[3,]    1    1    2    2
#[4,]    2    2    2    3

由于用于特征分解的算法不同，上述R因子与返回的因子不同。cholcovR 使用 LAPACK 例程DSYVER，其中进行了一些旋转以使特征值不增加。MATLABcholcov不是开源的，所以我不确定它使用什么算法。但是很容易证明它没有按非递增顺序排列特征值。

考虑T返回的因子cholcov：

T <- structure(c(-0.2113, 0.7887, 1.1547, 0.7887, -0.2113, 1.1547, 
-0.5774, -0.5774, 1.1547, 0, 0, 1.7321), .Dim = 3:4)

我们可以通过

rowSums(T ^ 2)
# [1] 1.000086 1.000086 7.000167

由于不精确，存在一些舍入误差T，但我们可以清楚地看到特征值是1, 1, 7。另一方面，我们有7, 1, 1来自 R (recall D1)。

于 2016-11-23T17:36:17.030 回答