Rの `chol`はMATLABの` cholcov`とは異なります。コレスキーに似た共分散分解を行うには？

Question

私はcholcov()を使ってMatlabで行われているように、Rでコレスキー的な共分散分解を再現しようとしています。例はhttps://uk.mathworks.com/help/stats/cholcov.htmlから取られます。彼らの例のように、元のcholcov()機能の

結果：私はR.でこのTを複製しようとしています

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

は、私が試した：

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のソリューションとはまったく異なります。私は、多分それは共分散行列のコレスキー組成だろうと思った：

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 

を取得します。

Matlabのcholcov()が計算されているので、どのようにしてRで複製できますか？

Answer 1

この場合、R関数cholを本質的に悪用しています。 MATLABのcholcov関数は複合関数です。

• 共分散が正の場合、コレスキー分解が行われ、フルランクの上三角コレスキー係数が返されます。
• 共分散が正の半正定値の場合は、固有分解を行い、長方形の行列を返します。

一方、Rからのcholは、Choleksy分解のみを行います。あなたが与えた例は、C1で、後者の場合に該当します。したがって

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

率ます：だから、我々は数値ランクが3であるので

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 

は、我々は最後の固有値と固有ベクトルをドロップR.にeigen機能に頼る必要があります欲しいものは：

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因子が原因固有分解のために使用される様々なアルゴリズムにcholcovによって返されるものと同じではありません。 RはLAPACKルーチンDSYVERを使用しており、固有値が増加しないようにピボットが行われます。MATLABのcholcovはオープ​​ンソースではないので、どのアルゴリズムを使用しているのかよくわかりません。しかし、それが非増加の順序で固有値を配列しないことを実証することは容易である。そこTが正確ではないので、いくつかの丸め誤差がありますが、我々ははっきりと見ることができる我々は

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

により固有値を得ることができます

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) 

：

はcholcovによって返さ要因Tを考えてみましょう

その固有値は 1, 1, 7です。一方、Rからの 7, 1, 1（リコール D1）があります。