R's `chol` отличается от MATLAB`cholcov`. Как сделать разложение ковариантности по Холецкому?

Question

Я пытался воспроизвести разложение ковариантности в R по типу холески - как это делается в 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 в Р. я попробовал:

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

что тоже не правильно.

Кто-нибудь может дать мне подсказку, как cholcov() в Matlab рассчитывается так, чтобы я мог повторить его в R?

Answer 1

Вы по сути злоупотребляете функцией R chol в этом случае. cholcov Функция от MATLAB является составной функцией.

• Если ковариация положительна, она выполняет факторизацию Холецкого, возвращая полный ранг верхнего треугольного фактора Холецкого;
• Если ковариация является положительно-полуопределенной, она выполняет собственное разложение, возвращая прямоугольную матрицу.

С другой стороны, chol из р только холексовая факторизация. Пример, который вы даете, C1, попадает во второй случай. Итак, мы должны прибегнуть к eigen функция в R.

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 коэффициент выше не такой, как тот, который возвращается cholcov из-за различных алгоритмов, используемых для собственной факторизации. R использует подпрограмму LAPACK DSYVER в котором некоторое вращение выполняется так, чтобы собственные значения не увеличивались. от Matlab cholcov не с открытым исходным кодом, поэтому я не уверен, какой алгоритм он использует. Но легко показать, что он не упорядочивает собственные значения в не возрастающем порядке.

Учитывайте фактор 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 (вспомните D1).