¿Cómo puedo generar datos con una matriz de correlación preespecificada?

19

Estoy tratando de generar una secuencia aleatoria correlacionada con media = 0 , varianza = , coeficiente de correlación = . En el siguiente código, uso & como desviaciones estándar y & como medio.0.810.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

Esto me da el correcto corrcoef()de 0.8 entre xy y. Mi pregunta es cómo puedo generar una serie de medios si quiero zque también esté correlacionado con y(con la misma correlación ), pero no con . ¿Hay una fórmula particular que necesito saber? Encontré uno pero no pude entenderlo.r=0.8x

anisa
fuente

Respuestas:

21

Parece que está preguntando cómo generar datos con una matriz de correlación particular.

Un hecho útil es que si tiene un vector aleatorio con matriz de covarianza Σ , entonces el vector aleatorio A x tiene media A E ( x ) y matriz de covarianza Ω = A Σ A T . Por lo tanto, si comienza con datos que tienen una media de cero, la multiplicación por A no cambiará eso, por lo que su primer requisito se cumple fácilmente. xΣAxAE(x)Ω=AΣATA

Digamos que usted comienza con los datos correlacionados (media cero) (es decir, la matriz de covarianza es diagonal) - ya que estamos hablando de la matriz de correlación, vamos a tomar . Puede transformar esto en datos con una matriz de covarianza dada eligiendo A como la raíz cuadrada cholesky de Ω ; entonces A x tendría la matriz de covarianza deseada Ω .Σ=IAΩAxΩ

En su ejemplo, parece querer algo como esto:

Ω=(1.80.81.80.81)

Desafortunadamente, esa matriz no es positiva definida, por lo que no puede ser una matriz de covarianza; puede verificar esto al ver que el determinante es negativo. Quizás, en cambio

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

bastaría. No estoy seguro de cómo calcular la raíz cuadrada de Cholesky en Matlab (que parece ser lo que estás usando), pero Rpuedes usar la chol()función.

En este ejemplo, para los dos s enumerados anteriormente, los múltiplos de matriz adecuados (respectivamente) seríanΩ

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

El Rcódigo utilizado para llegar a esto fue:

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
fuente
1
La función MATLAB también se llama chol. Tenga en cuenta que esto puede ser bastante inestable numéricamente si es casi singular. En ese caso, usar la raíz cuadrada simétrica obtenida, por ejemplo, a través de la SVD, es a menudo una mejor opción en términos de estabilidad numérica. :)Ω
cardenal
1
Por supuesto, eso es correcto @cardinal: muchas cosas teóricamente justificadas salen mal cuando intentas hacer cosas numéricamente con matrices casi singulares. Estaba (convenientemente) imaginando la situación en la que la matriz de correlación objetivo no estaba en el ámbito donde esto era un problema. Es bueno que hayas señalado esto - gracias (y gracias por la edición a mi otra respuesta)
Macro
1
La razón principal por la que estaba pensando en esto se debió a su buen ojo al reconocer que la primera sugerencia del OP ni siquiera fue definitiva y positiva. Y, con suerte, la edición de la otra pregunta no fue demasiado celosa; Me gustan ambas respuestas.
cardenal
7

Si está utilizando R, también puede usar la función mvrnorm del paquete MASS, suponiendo que desea variables normalmente distribuidas. La implementación es similar a la descripción de Macro anterior, pero utiliza los vectores propios de la matriz de correlación en lugar de la descomposición colesky y la escala con una descomposición de valor singular (si la opción empírica se establece en verdadero).

XΣ is a positive definite correlation matrix with eigenvectors γ, and λ is a square matrix with the square root eigen values from Σ along the diagonal then:

X=γλXT

Where X' is a normally distributed matrix with correlation matrix of Σ and column means are the same as X.

Note that the correlation matrix have to be positive definite, but converting it with the nearPD function from the Matrix package in R will be useful.

zzk
fuente
1

An alternative solution without cholesky factorization is the following. Let Σy the desired covariance matrix and suppose you have data x with Σx=I. Suppose Σy is positive definite with Λ the diagonal matrix of the eigenvalues and V the matrix of column eigenvectors .

You can write Σy=VΛVT=(VΛ)(ΛTVT)=AAT.

y=Ax generate the desired data.

Mario Sansone
fuente