¿Cómo generar números aleatorios correlacionados (dado los medios, varianzas y grado de correlación)?

Question

Lo siento si esto parece un poco básico, pero supongo que solo estoy buscando confirmar la comprensión aquí. Tengo la sensación de que tendría que hacer esto en dos pasos, y he comenzado a intentar entender las matrices de correlación, pero parece ser realmente complejo. Estoy buscando una explicación concisa (idealmente con sugerencias hacia una solución de seudocódigo) de una buena manera, idealmente rápida, de generar números aleatorios correlacionados.

Dados dos variables pseudoaleatorias altura y peso con medias y varianzas conocidas, y una correlación dada, creo que básicamente estoy tratando de entender cómo debería verse este segundo paso:

   altura = gaussianPdf(altura.media, altura.varianza)
   peso = gaussianPdf(media_correlacionada(altura.media, coeficiente_de_correlacion), 
                        varianza_correlacionada(altura.varianza, 
                        coeficiente_de_correlacion))

• ¿Cómo calculo la media y varianza correlacionadas? Pero quiero confirmar que ese es realmente el problema relevante aquí.
• ¿Necesito recurrir a la manipulación de matrices? ¿O tengo algo más muy erróneo en mi enfoque básico de este problema?

Answer 1

Para responder a tu pregunta sobre "una forma buena, idealmente rápida de generar números aleatorios correlacionados": Dada una matriz de varianza-covarianza deseada $C$ que, por definición, es definida positiva, la descomposición de Cholesky de la misma es: $C$=$LL^T$; siendo $L$ una matriz triangular inferior.

Si ahora usas esta matriz $L$ para proyectar un vector de variables aleatorias no correlacionadas $X$, la proyección resultante $Y = LX$ será la de variables aleatorias correlacionadas.

Puedes encontrar una explicación concisa de por qué esto sucede aquí.

Answer 2

+1 a @user11852 y, @jem77bfp, estas son buenas respuestas. Permíteme abordar esto desde una perspectiva diferente, no porque crea que necesariamente sea mejor en la práctica, sino porque creo que es instructivo. Aquí hay algunos hechos relevantes que ya conocemos:

1. $r$ es la pendiente de la línea de regresión cuando tanto $X$ como $Y$ están estandarizados, es decir, $\mathcal N(0,1)$,

2. $r^2$ es la proporción de la varianza en $Y$ atribuible a la varianza en $X,

   (también, de las reglas para varianzas):

3. la varianza de una variable aleatoria multiplicada por una constante es la constante al cuadrado veces la varianza original:
   $$\text{Var}[aX]=a^2\text{Var}[X]$$

4. las varianzas se suman, es decir, la varianza de la suma de dos variables aleatorias (asumiendo que son independientes) es la suma de las dos varianzas:
   $$\text{Var}[X+\varepsilon]=\text{Var}[X]+\text{Var}[\varepsilon]$$

Ahora, podemos combinar estos cuatro hechos para crear dos variables normales estándar cuyas poblaciones tendrán una correlación dada, $r$ (más correctamente, $\rho$), aunque las muestras que generes tendrán correlaciones de muestra que varían. La idea es crear una variable seudorandom, $X$, que sea normal estándar, $\mathcal N(0,1)$, y luego encontrar un coeficiente, $a$, y una varianza de error, $v_e$, de modo que $Y \sim\mathcal N(0,a^2+v_e)$, donde $a^2+v_e=1$. (Nota que $|a|$ debe ser $\le 1$ para que esto funcione, y además, $a=r$). Así que comienzas con el $r$ que deseas; ese es tu coeficiente, $a$. Luego descubres la varianza de error que necesitarás, que es $1-r^2$. (Si tu software te pide que uses la desviación estándar, toma la raíz cuadrada de ese valor). Finalmente, para cada variable seudorandom, $x_i$, que hayas generado, genera una variable de error seudorandom, $e_i$, con la varianza de error apropiada $v_e$, y calcula la variable seudorandom correlacionada, $y_i$, multiplicando y sumando.

Si deseas hacer esto en R, el siguiente código podría funcionar para ti:

correlatedValue = function(x, r){
  r2 = r**2
  ve = 1-r2
  SD = sqrt(ve)
  e  = rnorm(length(x), mean=0, sd=SD)
  y  = r*x + e
  return(y)
}

set.seed(5)
x = rnorm(10000)
y = correlatedValue(x=x, r=.5)

cor(x,y)
[1] 0.4945964

(Edit: Olvidé mencionar:) Como lo he descrito, este procedimiento te da dos variables correlacionadas normales estándar. Si no deseas normales estándar, sino que deseas que las variables tengan ciertos promedios específicos (no 0) y desviaciones estándar (no 1), puedes transformarlas sin afectar la correlación. Por lo tanto, restarías el promedio observado para asegurarte de que el promedio sea exactamente $0$, multiplicarías la variable por la desviación estándar que deseas y luego sumarías el promedio deseado. Si deseas que el promedio observado fluctúe normalmente alrededor del promedio deseado, agregarías la diferencia inicial nuevamente. Básicamente, esta es una transformación de z-score al revés. Debido a que esta es una transformación lineal, la variable transformada tendrá la misma correlación con la otra variable que antes.

Nuevamente, esto, en su forma más simple, solo te permite generar un par de variables correlacionadas (esto podría escalarse, pero se vuelve complicado rápidamente), y ciertamente no es la forma más conveniente de hacer el trabajo. En R, querrías usar ?mvrnorm en el paquete MASS, tanto porque es más fácil como porque puedes generar muchas variables con una matriz de correlación poblacional dada. Sin embargo, creo que vale la pena haber pasado por este proceso para ver cómo algunos principios básicos se desarrollan de una manera sencilla.

Answer 3

En general, esto no es algo simple de hacer, pero creo que hay paquetes para la generación de variables normales multivariadas (al menos en R, ver mvrnorm en el paquete MASS), donde solo tienes que ingresar una matriz de covarianza y un vector de medias.

También hay otro enfoque "constructivo". Supongamos que queremos modelar un vector aleatorio $(X_1,X_2)$ y tenemos su función de distribución $F(x_1,x_2)$. El primer paso es obtener la función de distribución marginal; es decir, integrar $F$ sobre todos los $x_2$: $$F_{X_1}(x_1)= \int_{-\infty}^{\infty} F(x_1,x_2)dx_2.$$ Luego encontramos $F^{-1}_{X_1}$ - la función inversa de $F_{X_1}$ - e insertamos una variable aleatoria $\xi_1$ que está uniformemente distribuida en el intervalo $[0,1]$. En este paso generamos la primera coordenada $\hat{x}_1=F^{-1}_{X_1}(\xi)$.

Ahora, una vez que tenemos una coordenada, necesitamos insertarla en nuestra función de distribución inicial $F(x_1,x_2)$ y luego obtener una función de distribución condicional con la condición $x_1=\hat{x}_1$: $$F(x_2 | X_1=\hat{x}_1)= \frac{F(\hat{x}_1,x_2)}{f_{X_1}(\hat{x}_1)},$$ donde $f_{X_1}$ es una función de densidad de probabilidad de la distribución marginal $X_1$; es decir, $F'_{X_1}(x_1)=f_{X_1}(x_1)$.

Luego nuevamente generas una variable distribuida de forma uniforme $\xi_2$ en $[0,1]$ (independiente de $\xi_1$) e insertas en la inversa de $F(x_2 | X_1=\hat{x}_1)$. Por lo tanto, obtienes $\hat{x}_2=(F(x_2 | X_1=\hat{x}_1))^{-1}(\xi)$; es decir, $\hat x_2$ satisface $F(\hat x_2 | X_1=\hat{x}_1) = \xi$. Este método se puede generalizar a vectores con más dimensiones, pero su desventaja es que tienes que calcular, analítica o numéricamente, muchas funciones. La idea también se puede encontrar en este artículo: http://www.econ-pol.unisi.it/dmq/pdf/DMQ_WP_34.pdf.

Si no entiendes el significado de insertar una variable uniforme en una función de distribución de probabilidad inversa, intenta hacer un esquema del caso univariante y luego recuerda cuál es la interpretación geométrica de la función inversa.

Answer 4

Usando la descomposición de Cholesky, creo que esta función podría funcionar en R:

Cholesky_samples <- function(r, n, mean1, mean2, sd1, sd2){
C <- matrix(c(1,r,r,1), nrow = 2)
L <- chol(C)
X1 <- rnorm(n)
X2 <- rnorm(n)
X <- rbind(X1,X2)

Y <- t(L)%*%X
Y1 <- Y[1,]
Y2 <- Y[2,]

N_1 <<- Y[1,] * sd1 + mean1
N_2 <<- Y[2,] * sd2 + mean2

sample_measures <<- as.data.frame(c(mean(N_1), mean(N_2), sd(N_1), sd(N_2), cor(N_1, N_2)), 
                  names<-c("mean N_1", "mean N_2", "SD N_1", "SD N_2","cor(N_1,N_2)"))
sample_measures
}

Sin embargo, de manera más flexible en cuanto a las distribuciones marginales de las variables correlacionadas devueltas, creo que podemos usar la transformación Uniforme, la cual he integrado en la siguiente función:

cor_samples <- function(r, n, mean1, mean2, sd1, sd2){
  C <- matrix(c(1,r,r,1), nrow = 2)
  require(mvtnorm)
  SN <- rmvnorm(mean = c(0,0), sig = C, n = n)
  U <- pnorm(SN)
  U1 <- U[,1]
  U2 <- U[,2]

  Y1 <<- qnorm(U1, mean = mean1,sd = sd1) 
  Y2 <<- qnorm(U2, mean = mean2,sd = sd2) 

  sample_measures <<- as.data.frame(c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1,Y2)), 
                      names<-c("mean Y1", "mean Y2", "SD Y1", "SD Y2", "cor(Y1,Y2)"))
  sample_measures
}

En el primer método, la salida incluye dos muestras de un tamaño seleccionado, correlacionadas como se prescribe, y distribuidas normalmente.

En contraste, el segundo método es fácil de adaptar para producir muestras correlacionadas con distribuciones marginales seleccionadas normales o no normales. Por ejemplo, podríamos obtener dos muestras correlacionadas siguiendo Y1 ~ $Beta(\alpha=\beta=0.5)$ y Y2 ~ $Exp(\lambda=0.5)$ todavía con una aproximación razonable a la correlación elegida (digamos, r = 0.7):

cor_samples <- function(r, n){
  C <- matrix(c(1,r,r,1), nrow = 2)
  require(mvtnorm)
  SN <- rmvnorm(mean = c(0,0), sig = C, n = n)
  U <- pnorm(SN)
  U1 <- U[,1]
  U2 <- U[,2]

  Y1 <<- qbeta(U1, 0.5,0.5) 
  Y2 <<- qexp(U2, 0.5) 

  sample_measures <<- as.data.frame(c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1,Y2)), 
                      names<-c("mean Y1", "mean Y2", "SD Y1", "SD Y2", "cor(Y1,Y2)"))
  sample_measures
}

Los histogramas de la salida serían:

Answer 5

Si estás listo para renunciar a la eficiencia, puedes usar un algoritmo desechable. Su ventaja es que permite cualquier tipo de distribuciones (no solo la Gaussiana).

Comienza generando dos secuencias no correlacionadas de números aleatorios $\{x_i\}_{i=1}^N$ y $\{y_i\}_{i=1}^N$ con las distribuciones deseadas. Sea $C$ el valor deseado del coeficiente de correlación. Luego haz lo siguiente:

1) Calcula el coeficiente de correlación $c_{old}=corr(\{x_i\},\{y_i\})$

2) Genera dos números aleatorios $n_1$ y $n_2: 1 \leq n_{1,2} \leq N$

3) Intercambia los números $x_{n_1}$ y $x_{n_2}$

4) Calcula la nueva correlación $c_{new}=corr( \{x_i\},\{y_i\})$

5) Si $|C-c_{new}| < |C-c_{old}|$ entonces mantén el intercambio. De lo contrario, deshaz el intercambio.

6) Si $|C-c| < \epsilon$ entonces detente, de lo contrario regresa al paso 1)

Los intercambios aleatorios no alterarán la distribución marginal de ${x_i}$.

¡Buena suerte!