Как создать произвольную ковариационную матрицу

21

Например, в R, MASS::mvrnorm()функция полезна для генерации данных, чтобы продемонстрировать различные вещи в статистике. Он принимает обязательный Sigmaаргумент, который представляет собой симметричную матрицу, определяющую ковариационную матрицу переменных. Как бы я создал симметричную матрицу с произвольными записями?n×n

RSL
источник
3
Я думаю, что этот вопрос выиграл бы от редактирования, чтобы сосредоточиться на том, «как я могу создать произвольную ковариационную матрицу», а не на аспекте кодирования. Здесь, безусловно, есть основополагающая статистическая проблема, о чем свидетельствует ответ.
Серебряная рыба

Ответы:

22

Создать матрицу с произвольными значениямиn×nA

и затем используйте качестве вашей ковариационной матрицы. Σ=ATA

Например

n <- 4  
A <- matrix(runif(n^2)*2-1, ncol=n) 
Sigma <- t(A) %*% A
Генри
источник
Точно так же Sigma <- A + t(A).
RSL
6
@MoazzemHossen: Ваше предложение создаст симметричную матрицу, но оно не всегда может быть положительным полуопределенным (например, ваше предложение может привести к матрице с отрицательными собственными значениями), и поэтому оно может не подходить в качестве ковариационной матрицы
Генри
Да, я заметил, что R возвращает ошибку в том случае, если мой предложенный способ дал неподходящую матрицу.
RSL
4
Обратите внимание, что если вы предпочитаете корреляционную матрицу для лучшей интерпретации, есть функция ? Cov2cor , которая может быть применена впоследствии.
gung - Восстановить Монику
1
@ B11b: Ваша ковариационная матрица должна быть положительной полуопределенной. Это наложило бы некоторые ограничения на значения ковариации, не совсем очевидные, когда n>2
Генри
24

Мне нравится контролировать объекты, которые я создаю, даже когда они могут быть произвольными.

Тогда рассмотрим, что все возможные ковариационные матрицы Σ можно выразить в видеn×nΣ

Σ=P Diagonal(σ1,σ2,,σn) P

где - ортогональная матрица и σ 1σ 2σ n0 .Pσ1σ2σn0

Геометрически это описывает ковариационную структуру с диапазоном главных компонент размеров . Эти компоненты указывают в направлениях строк P . См. Рисунки в разделе « Осмысление анализа главных компонент, собственных векторов и собственных значений» для примеров с n = 3 . Установка σ i будет устанавливать величины ковариаций и их относительные размеры, тем самым определяя любую желаемую форму эллипсоида. Ряды P ориентируют оси формы так, как вы предпочитаете.σiPn=3σiP

Одно алгебраическое и вычислительное преимущество этого подхода состоит в том, что когда , Σ легко инвертируется (что является обычной операцией на ковариационных матрицах):σn>0Σ

Σ1=P Diagonal(1/σ1,1/σ2,,1/σn) P.

σin2n

n <- 5
p <- qr.Q(qr(matrix(rnorm(n^2), n)))

n

ΣPσicrossprodRσ=(σ1,,σ5)=(5,4,3,2,1)

Sigma <- crossprod(p, p*(5:1))

σP

svd(Sigma)

Sigmaσ

Tau <- crossprod(p, p/(5:1))

zapsmall(Sigma %*% Tau)n×nσi01/σiσi

Whuber
источник
P
1
Стоит упомянуть, что единственные значения в svd(Sigma)будут переупорядочены - это смутило меня на минуту.
FrankD
1

Вы можете моделировать случайные положительно определенные матрицы из дистрибутива Wishart, используя функцию "rWishart" из широко используемого пакета "stats".

n <- 4
rWishart(1,n,diag(n))
Карлос Льоса
источник
1

Для этого есть специальная посылка clusterGeneration(написанная, в частности, Гарри Джо, известным в этой области).

Есть две основные функции:

  • genPositiveDefMat сгенерировать ковариационную матрицу, 4 разных метода
  • rcorrmatrix : создать матрицу корреляции

Быстрый пример:

library(clusterGeneration)
#> Loading required package: MASS
genPositiveDefMat("unifcorrmat",dim=3)
#> $egvalues
#> [1] 15.408962  5.673916  1.228842
#> 
#> $Sigma
#>          [,1]     [,2]     [,3]
#> [1,] 6.714871 1.643449 6.530493
#> [2,] 1.643449 6.568033 2.312455
#> [3,] 6.530493 2.312455 9.028815
genPositiveDefMat("eigen",dim=3)
#> $egvalues
#> [1] 8.409136 4.076442 2.256715
#> 
#> $Sigma
#>            [,1]       [,2]      [,3]
#> [1,]  2.3217300 -0.1467812 0.5220522
#> [2,] -0.1467812  4.1126757 0.5049819
#> [3,]  0.5220522  0.5049819 8.3078880

Создано в 2019-10-27 пакетом представлением (v0.3.0)

Наконец, обратите внимание, что альтернативный подход состоит в том, чтобы сначала попробовать с нуля, а затем использовать, Matrix::nearPD()чтобы сделать вашу матрицу положительно определенной.

Matifou
источник