Comment générer des données avec une matrice de corrélation prédéfinie?

19

J'essaie de générer une séquence aléatoire corrélée avec une moyenne = , une variance = , un coefficient de corrélation = . Dans le code ci-dessous, j'utilise & comme écart-type et & comme moyen.1 0,8010,8s1s2m1m2

p = 0.8 
u = randn(1, n)
v = randn(1, n)
x = s1 * u + m1
y = s2 * (p * u + sqrt(1 - p^2) * v) + m2

Cela me donne le correct corrcoef()de 0,8 entre xet y. Ma question est de savoir comment puis-je générer une série signifie si je veux zque cela soit également corrélé avec y(avec la même corrélation ), mais pas avec . Existe-t-il une formule particulière que je dois connaître? J'en ai trouvé un mais je ne pouvais pas le comprendre.r=0,8x

anisa
la source

Réponses:

21

Il semble que vous vous demandiez comment générer des données avec une matrice de corrélation particulière.

Un fait est utile que si vous avez un vecteur aléatoire avec la matrice de covariance Σ , le vecteur aléatoire A x a moyenne A E ( x ) et covariance matrice Ω = A Σ A T . Donc, si vous commencez avec des données qui ont une moyenne de zéro, la multiplication par A ne changera pas cela, donc votre première exigence est facilement satisfaite. XΣUNEXUNEE(X)Ω=UNEΣUNETUNE

Disons que vous commencez avec (zéro) moyenne des données non corrélées (la matrice de covariance est diagonale) - puisque nous parlons de la matrice de corrélation, nous allons juste prendre . Vous pouvez transformer cela en données avec une matrice de covariance donnée en choisissant A pour être la racine carrée cholesky de Ω - alors A x aurait la matrice de covariance Ω souhaitée .Σ=jeUNEΩUNEXΩ

Dans votre exemple, vous semblez vouloir quelque chose comme ceci:

Ω=(1.80.81.80.81)

Malheureusement, cette matrice n'est pas définie positive, elle ne peut donc pas être une matrice de covariance - vous pouvez le vérifier en voyant que le déterminant est négatif. Peut-être, à la place

Ω=(1.8.3.81.8.3.81)    or   Ω=(12/302/312/302/31)

suffirait. Je ne sais pas comment calculer la racine carrée cholesky dans matlab (qui semble être ce que vous utilisez) mais Rvous pouvez utiliser la chol()fonction.

Dans cet exemple, pour les deux énumérés ci-dessus, les multiples de matrice appropriés (respectivement) seraientΩ

UNE=(100.8.60.3.933.1972)    or   UNE=(1002/3.745300.8944.4472)

Le Rcode utilisé pour y parvenir était:

x = matrix(0,3,3)
x[1,]=c(1,.8,.3)
x[2,]=c(.8,1,.8)
x[3,]=c(.3,.8,1)
t(chol(x))

     [,1]      [,2]      [,3]
[1,]  1.0 0.0000000 0.0000000
[2,]  0.8 0.6000000 0.0000000
[3,]  0.3 0.9333333 0.1972027

x[1,]=c(1,2/3,0)
x[2,]=c(2/3,1,2/3)
x[3,]=c(0,2/3,1)
t(chol(x))

      [,1]      [,2]      [,3]
[1,] 1.0000000 0.0000000 0.0000000
[2,] 0.6666667 0.7453560 0.0000000
[3,] 0.0000000 0.8944272 0.4472136
Macro
la source
1
cholΩ
1
Bien sûr, c'est vrai @cardinal - beaucoup de choses théoriquement justifiées vont mal lorsque vous essayez de faire des choses numériquement avec des matrices presque singulières. J'imaginais (commodément) la situation où la matrice de corrélation cible n'était pas dans le domaine où c'était un problème. C'est bien que vous l'ayez signalé - merci (et merci pour la modification de mon autre réponse)
Macro
1
La principale raison pour laquelle j'y pensais était due à votre œil attentif pour reconnaître que la première suggestion du PO n'était même pas définie de manière positive. Et, espérons-le, la modification de l'autre question n'était pas trop zélée; J'aime ces deux réponses.
cardinal
7

Si vous utilisez R, vous pouvez également utiliser la fonction mvrnorm du package MASS, en supposant que vous souhaitiez des variables normalement distribuées. L'implémentation est similaire à la description de Macro ci-dessus, mais utilise les vecteurs propres de la matrice de corrélation au lieu de la décomposition et de la mise à l'échelle cholesky avec une décomposition en valeurs singulières (si l'option empirique est définie sur true).

XΣγλΣ

X=γλXT

ΣX

Notez que la matrice de corrélation doit être définie positive, mais la convertir avec la fonction nearPD du package Matrix dans R sera utile.

zzk
la source
1

ΣyXΣX=jeΣyΛV

Σy=VΛVT=(VΛ)(ΛTVT)=UNEUNET

y=UNEX

Mario Sansone
la source