Generar una variable aleatoria con una correlación definida con una variable o variables existentes

Para un estudio de simulación tengo para generar variables aleatorias que muestran un (población) de correlación prefined a una variable existente .$Y$

Miré en los Rpaquetes copulay CDVineque pueden producir distribuciones aleatorias multivariadas con una estructura de dependencia dada. Sin embargo, no es posible arreglar una de las variables resultantes a una variable existente.

Cualquier idea y enlaces a las funciones existentes son apreciados!

Conclusión: surgieron dos respuestas válidas, con diferentes soluciones:

1. Un R guión de caracal, que calcula una variable aleatoria con una correlación exacta (muestra) con una variable predefinida
2. Una R función que encontré yo mismo, que calcula una variable aleatoria con una correlación de población definida con una variable predefinida

[Adición de @ttnphns: me tomé la libertad de ampliar el título de la pregunta de un caso de variable fija única a un número arbitrario de variables fijas; es decir, cómo generar una variable que tenga correcciones predefinidas con algunas variables fijas y existentes]

Aquí hay otro: para los vectores con media 0, su correlación es igual al coseno de su ángulo. Entonces, una forma de encontrar un vector con exactamente la correlación deseada , correspondiente a un ángulo :r $x$$r$$\theta$

1. obtener el vector fijo y un vector aleatoriox $x_1$$x_2$
2. centrar ambos vectores (media 0), dando vectores , ˙ x$\dot{x}_{1}$$\dot{x}_{2}$
3. hacer ortogonal a (proyección sobre el subespacio ortogonal), dando ˙ x 1 ˙ x ⊥ $\dot{x}_{2}$$\dot{x}_{1}$$\dot{x}_{2}^{\perp}$
4. escala y a la longitud 1, dando y ˙ x ⊥ 2 ˉ x 1 ˉ x ⊥ $\dot{x}_{1}$$\dot{x}_{2}^{\perp}$$\bar{x}_{1}$$\bar{x}_{2}^{\perp}$
5. ˉ x 1θ ˉ x 1rx$\bar{x}_{2}^{\perp} + (1/\tan(\theta)) \cdot \bar{x}_{1}$ es el vector cuyo ángulo a es , y cuya correlación con es así . Esta es también la correlación con ya que las transformaciones lineales dejan la correlación sin cambios.$\bar{x}_{1}$$\theta$$\bar{x}_{1}$$r$$x_1$

Aquí está el código:

n     <- 20                    # length of vector
rho   <- 0.6                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- rnorm(n, 1, 1)        # fixed given data
x2    <- rnorm(n, 2, 0.5)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho

Para la proyección ortogonal , utilicé la descomposición para mejorar la estabilidad numérica, desde entonces simplemente .Q R P = Q Q $P$$QR$$P = Q Q'$

Describiré la solución más general posible. Resolver el problema en esta generalidad nos permite lograr una implementación de software notablemente compacta: solo Rbastan dos líneas cortas de código.

Elija un vector , de la misma longitud que , de acuerdo con la distribución que desee. Deje que ser los residuos de la regresión de mínimos cuadrados de contra : este extrae el componente de . Mediante la adición de nuevo un múltiplo adecuado de a , podemos producir un vector que tiene cualquier correlación deseada con . Hasta una constante aditiva arbitraria y una constante multiplicativa positiva, que puede elegir de cualquier manera, la solución esY Y ⊥ X Y Y X Y Y ⊥ ρ $X$$Y$$Y^\perp$$X$$Y$$Y$$X$$Y$$Y^\perp$$\rho$$Y$

(" " representa cualquier cálculo proporcional a una desviación estándar).$\operatorname{SD}$

Aquí está el Rcódigo de trabajo . Si no proporciona , el código extraerá sus valores de la distribución Normal estándar multivariante.$X$

complement <- function(y, rho, x) {
  if (missing(x)) x <- rnorm(length(y)) # Optional: supply a default if `x` is not given
  y.perp <- residuals(lm(x ~ y))
  rho * sd(y.perp) * y + y.perp * sd(y) * sqrt(1 - rho^2)
}

Para ilustrar, me genera un aleatoria con componentes y produje que tiene diversas correlaciones especificado con esta . Todos fueron creados con el mismo vector inicial . Aquí están sus diagramas de dispersión. Los "rugplots" en la parte inferior de cada panel muestran el vector común .50 X Y ; ρ Y X = ( 1 , 2 , … , 50 ) $Y$$50$$X_{Y;\rho}$$Y$$X=(1,2,\ldots, 50)$$Y$

Hay una notable similitud entre las parcelas, no está allí :-).

Si desea experimentar, aquí está el código que produjo estos datos y la figura. (No me molesté en usar la libertad de cambiar y escalar los resultados, que son operaciones fáciles).

y <- rnorm(50, sd=10)
x <- 1:50 # Optional
rho <- seq(0, 1, length.out=6) * rep(c(-1,1), 3)
X <- data.frame(z=as.vector(sapply(rho, function(rho) complement(y, rho, x))),
                rho=ordered(rep(signif(rho, 2), each=length(y))),
                y=rep(y, length(rho)))

library(ggplot2)
ggplot(X, aes(y,z, group=rho)) + 
  geom_smooth(method="lm", color="Black") + 
  geom_rug(sides="b") + 
  geom_point(aes(fill=rho), alpha=1/2, shape=21) +
  facet_wrap(~ rho, scales="free")

Por cierto, este método se generaliza fácilmente a más de una : si es matemáticamente posible, encontrará una con correlaciones específicas conjunto de . Simplemente use los mínimos cuadrados ordinarios para eliminar los efectos de todos los de y forme una combinación lineal adecuada de y los residuos. (Ayuda hacer esto en términos de una base dual para , que se obtiene calculando un pseudo-inverso. El siguiente código usa la SVD de para lograr eso).X Y 1 , Y 2 , ... , Y k ; ρ 1 , ρ 2 , … , ρ k Y i Y i X Y i Y $Y$$X_{Y_1,Y_2,\ldots,Y_k;\rho_1,\rho_2,\ldots,\rho_k}$$Y_i$$Y_i$$X$$Y_i$$Y$$Y$

Aquí hay un boceto del algoritmo R, donde se dan como columnas de una matriz :$Y_i$y

y <- scale(y)             # Makes computations simpler
e <- residuals(lm(x ~ y)) # Take out the columns of matrix `y`
y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
return(y.dual %*% rho + sqrt(sigma2)*e)

La siguiente es una implementación más completa para aquellos que deseen experimentar.

complement <- function(y, rho, x) {
  #
  # Process the arguments.
  #
  if(!is.matrix(y)) y <- matrix(y, ncol=1)
  if (missing(x)) x <- rnorm(n)
  d <- ncol(y)
  n <- nrow(y)
  y <- scale(y) # Makes computations simpler
  #
  # Remove the effects of `y` on `x`.
  #
  e <- residuals(lm(x ~ y))
  #
  # Calculate the coefficient `sigma` of `e` so that the correlation of
  # `y` with the linear combination y.dual %*% rho + sigma*e is the desired
  # vector.
  #
  y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
  sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
  #
  # Return this linear combination.
  #
  if (sigma2 >= 0) {
    sigma <- sqrt(sigma2) 
    z <- y.dual %*% rho + sigma*e
  } else {
    warning("Correlations are impossible.")
    z <- rep(0, n)
  }
  return(z)
}
#
# Set up the problem.
#
d <- 3           # Number of given variables
n <- 50          # Dimension of all vectors
x <- 1:n         # Optionally: specify `x` or draw from any distribution
y <- matrix(rnorm(d*n), ncol=d) # Create `d` original variables in any way
rho <- c(0.5, -0.5, 0)          # Specify the correlations
#
# Verify the results.
#
z <- complement(y, rho, x)
cbind('Actual correlations' = cor(cbind(z, y))[1,-1],
      'Target correlations' = rho)
#
# Display them.
#
colnames(y) <- paste0("y.", 1:d)
colnames(z) <- "z"
pairs(cbind(z, y))

Aquí hay otro enfoque computacional (la solución está adaptada de una publicación del foro de Enrico Schumann). Según Wolfgang (ver comentarios), esto es computacionalmente idéntico a la solución propuesta por ttnphns.

A diferencia de la solución de caracal, no produce una muestra con la correlación exacta de , sino dos vectores cuya correlación de población es igual a .$\rho$$\rho$

La siguiente función puede calcular una distribución de muestra bivariada extraída de una población con un determinado . Calcula dos variables aleatorias o toma una variable existente (pasada como parámetro ) y crea una segunda variable con la correlación deseada:$\rho$x

# returns a data frame of two variables which correlate with a population correlation of rho
# If desired, one of both variables can be fixed to an existing variable by specifying x
getBiCop <- function(n, rho, mar.fun=rnorm, x = NULL, ...) {
     if (!is.null(x)) {X1 <- x} else {X1 <- mar.fun(n, ...)}
     if (!is.null(x) & length(x) != n) warning("Variable x does not have the same length as n!")

     C <- matrix(rho, nrow = 2, ncol = 2)
     diag(C) <- 1

     C <- chol(C)

     X2 <- mar.fun(n)
     X <- cbind(X1,X2)

     # induce correlation (does not change X1)
     df <- X %*% C

     ## if desired: check results
     #all.equal(X1,X[,1])
     #cor(X)

     return(df)
}

La función también puede usar distribuciones marginales no normales ajustando el parámetro mar.fun. Tenga en cuenta, sin embargo, que arreglar una variable solo parece funcionar con una variable normalmente distribuida x. (que podría estar relacionado con el comentario de Macro).

También tenga en cuenta que el "pequeño factor de corrección" de la publicación original se eliminó, ya que parece sesgar las correlaciones resultantes, al menos en el caso de las distribuciones gaussianas y las correlaciones de Pearson (también ver comentarios).

Deje que sea ​​su variable fija y desee generar una variable que se correlacione con en la cantidad . Si está estandarizado, entonces (porque es el coeficiente beta en regresión simple) , donde es una variable aleatoria de la distribución normal que tiene media y . La correlación observada entre los datos e será aproximadamente ; e pueden verse como muestras aleatorias de población normal bivariada (si$X$$Y$$X$$r$$X$$r$$Y= rX+E$$E$$0$$\text{sd}=\sqrt{1-r^2}$$X$$Y$$r$$X$$Y$$X$ es de lo normal) con .$\rho=r$

Ahora bien, si se quiere alcanzar la correlación de dos variables en la muestra exactamente , es necesario disponer que tiene cero correlación con . Este ajuste a cero se puede alcanzar modificando iterativamente. Bueno, con solo dos variables, una dada ( ) y otra para generar ( ), el número suficiente de iteraciones es en realidad 1, pero con múltiples variables dadas ( ) se necesitarán iteraciones.$r$$E$$X$$E$$X$$Y$$X_1, X_2, X_3,...$

Cabe señalar que si es normal, en el primer procedimiento ("aproximado ") también será normal; sin embargo, en el ajuste iterativo de a la " exacta " es probable que pierda normalidad porque el ajuste explota los valores de los casos de forma selectiva.$X$$r$$Y$$Y$$r$$Y$

Actualización 11 de noviembre de 2017. Hoy me he encontrado con este viejo hilo y decidí ampliar mi respuesta mostrando el algoritmo del ajuste iterativo sobre el que estaba hablando inicialmente.

Aquí hay una solución iterativa sobre cómo entrenar una variable aleatoriamente simulada o preexistente para correlacionar o covariar exactamente como lo deseamos (o muy cerca de eso, dependiendo del número de iteraciones) con un conjunto de variables dadas s (estas no pueden modificarse).$Y$ $X$

Descargo de responsabilidad: esta solución iterativa que he encontrado es inferior a la excelente basada en encontrar la base dual y propuesta por @whuber en este hilo hoy. La solución de @ whuber no es iterativa y, lo que es más importante para mí, parece estar afectando los valores de la variable de entrada "pig" algo menos que el algoritmo "my" (sería una ventaja si la tarea es "corregir" la variable existente y no generar una variable aleatoria desde cero). Aún así, estoy publicando el mío por curiosidad y porque funciona (véase también la nota al pie).

Entonces, hemos dado las variables (fijas) , y la variable que es un "cerdo" de valores generado aleatoriamente o es una variable de datos existente cuyos valores necesitamos "corregir" - para obtener exactamente a las correlaciones (o pueden ser covarianzas) con las s. Todos los datos deben ser continuos; en otras palabras, debe haber una gran cantidad de valores únicos.$X_1, X_2,...,X_m$$Y$$Y$$r_1, r_2,...,r_m$$X$

La idea: realizar un ajuste iterativo de residuos. Conociendo las correlaciones / covarianzas deseadas (objetivo), podemos calcular los valores pronosticados para la utilizando las s como predictores lineales múltiples. Después de obtener los residuos iniciales (de la actual y la predicción ideal), entrénelos iterativamente para que no se correlacionen con los predictores. Al final, recupere con los residuos. (El procedimiento fue mi propio invento experimental de la rueda hace muchos años cuando no conocía nada de la teoría; luego lo codifiqué en SPSS).$Y$$X$$Y$$Y$

1. Convierta el objetivo s en sumas de productos cruzados multiplicándolos por : . ( es un índice variable ).$r$$\text{df}=n-1$$S_j=r_j \text{df}$$j$$X$

2. Estandarizar en Z todas las variables (centrar cada una, luego dividir por la desviación st. Calculada en el anterior ). y s son, por lo tanto, estándar. Las sumas de cuadrados observadas son ahora = .$\text{df}$$Y$$X$$\text{df}$

3. Calcular los coeficientes de predicción de regressional por s de acuerdo con el objetivo s: .$Y$$X$$r$$\bf b=(X'X)^{-1} S$

4. Calcule los valores pronosticados para : .$Y$$\hat{Y}=\bf Xb$

5. Calcule los residuos .$E=Y-\hat{Y}$

6. Calcule la suma de cuadrados (objetivo) necesaria para los residuos: .$SS_S=\text{df}-SS_{\hat {Y}}$

7. (Comience a repetir.) Calcule las sumas observadas de productos cruzados entre actual y cada :$E$$X_j$$C_j= \sum_{i=1}^n E_i X_{ij}$

8. Corrija los valores de con el objetivo de acercar todos los s a ( es un índice de caso):$E$$C$$0$$i$

   $$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m C_j X_{ij}} {n\sum_{j=1}^m X_{ij}^2}$$

   (el denominador no cambia en las iteraciones, calcule por adelantado)

   O, alternativamente, una fórmula más eficiente asegura además que la media de convierte en . Primero, centre en cada cálculo previo de la s en el paso 7, luego en este paso 8 corrija como:$E$$0$ $E$$C$

   $$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m \frac{C_j X_{ij}^3}{\sum_{i=1}^n X_{ij}^2}} {\sum_{j=1}^m X_{ij}^2}$$

   (de nuevo, los denominadores se conocen de antemano)$^1$

9. Traiga a su valor objetivo:$SS_E$$E_i[\text{corrected}]=E_i \sqrt{SS_S/SS_E}$

   Vaya al paso 7. (Haga, digamos, 10-20 iteraciones; cuanto mayor sea más iteraciones podrían necesitarse. Si el objetivo s fuera realista, es positivo, y si el tamaño de la muestra no es demasiado pequeño, las iteraciones siempre directo a la convergencia. Fin de iteración.)$m$$r$$SS_S$$n$

10. Listo: Todos los s son casi cero ahora, lo que significa que los residuos han sido entrenados para restaurar objetivo . Calcular el ajuste : .E r Y Y [ corregido ] = Y + $C$$E$$r$$Y$$Y[\text{corrected}]=\hat{Y}+E$

11. La obtenida está casi estandarizada. Como último golpe, es posible que desee estandarizarlo con precisión, de nuevo como lo hizo en el paso 2.$Y$

12. Puede suministrar a cualquier variación y decir lo que quiera. En realidad, entre las cuatro estadísticas: min , max , mean , st. dev . - puede seleccionar cualquiera de los dos valores y transformar linealmente la variable para que los posea sin alterar las s (correlaciones) que ha alcanzado (todo se llama reescalamiento lineal).$Y$$r$

Para advertir nuevamente lo que se dijo anteriormente. Con ese tirón de exactamente a la , la salida no tiene que distribuirse normalmente.r $Y$$r$$Y$

Y $^1$ La fórmula de corrección puede ser aún más sofisticada, por ejemplo, para asegurar una mayor homocedasticidad (en términos de sumas de cuadrados) de con cada también, simultáneamente con el logro de las correlaciones, - he implementado un código para eso también. (No sé si esa tarea "doble" se puede resolver a través de un enfoque más ordenado, no narrativo, como el de Whuber ).$Y$$X$

Sentí ganas de programar un poco, así que tomé la respuesta eliminada de @ Adam y decidí escribir una buena implementación en R. Me concentro en usar un estilo funcionalmente orientado (es decir, un bucle de estilo rápido). La idea general es tomar dos vectores, permutar aleatoriamente uno de los vectores hasta que se haya alcanzado una cierta correlación entre ellos. Este enfoque es de fuerza bruta, pero es simple de implementar.

Primero creamos una función que permuta aleatoriamente el vector de entrada:

randomly_permute = function(vec) vec[sample.int(length(vec))]
randomly_permute(1:100)
  [1]  71  34   8  98   3  86  28  37   5  47  88  35  43 100  68  58  67  82
 [19]  13   9  61  10  94  29  81  63  14  48  76   6  78  91  74  69  18  12
 [37]   1  97  49  66  44  40  65  59  31  54  90  36  41  93  24  11  77  85
 [55]  32  79  84  15  89  45  53  22  17  16  92  55  83  42  96  72  21  95
 [73]  33  20  87  60  38   7   4  52  27   2  80  99  26  70  50  75  57  19
 [91]  73  62  23  25  64  51  30  46  56  39

... y crea algunos datos de ejemplo

vec1 = runif(100)
vec2 = runif(100)

... escriba una función que permute el vector de entrada y lo correlacione con un vector de referencia:

permute_and_correlate = function(vec, reference_vec) {
    perm_vec = randomly_permute(vec)
    cor_value = cor(perm_vec, reference_vec)
    return(list(vec = perm_vec, cor = cor_value))
  }
permute_and_correlate(vec2, vec1)
$vec
  [1] 0.79072381 0.23440845 0.35554970 0.95114398 0.77785348 0.74418811
  [7] 0.47871491 0.55981826 0.08801319 0.35698405 0.52140366 0.73996913
 [13] 0.67369873 0.85240338 0.57461506 0.14830718 0.40796732 0.67532970
 [19] 0.71901990 0.52031017 0.41357545 0.91780357 0.82437619 0.89799621
 [25] 0.07077250 0.12056045 0.46456652 0.21050067 0.30868672 0.55623242
 [31] 0.84776853 0.57217746 0.08626022 0.71740151 0.87959539 0.82931652
 [37] 0.93903143 0.74439384 0.25931398 0.99006038 0.08939812 0.69356590
 [43] 0.29254936 0.02674156 0.77182339 0.30047034 0.91790830 0.45862163
 [49] 0.27077191 0.74445997 0.34622648 0.58727094 0.92285322 0.83244284
 [55] 0.61397396 0.40616274 0.32203732 0.84003379 0.81109473 0.50573325
 [61] 0.86719899 0.45393971 0.19701975 0.63877904 0.11796154 0.26986325
 [67] 0.01581969 0.52571331 0.27087693 0.33821824 0.52590383 0.11261002
 [73] 0.89840404 0.82685046 0.83349287 0.46724807 0.15345334 0.60854785
 [79] 0.78854984 0.95770015 0.89193212 0.18885955 0.34303707 0.87332019
 [85] 0.08890968 0.22376395 0.02641979 0.43377516 0.58667068 0.22736077
 [91] 0.75948043 0.49734797 0.25235660 0.40125309 0.72147500 0.92423638
 [97] 0.27980561 0.71627101 0.07729027 0.05244047

$cor
[1] 0.1037542

... e iterar mil veces:

n_iterations = lapply(1:1000, function(x) permute_and_correlate(vec2, vec1))

Tenga en cuenta que reglas de alcance de R asegurar que vec1y vec2se encuentran en el medio ambiente mundial, fuera de la función anónima utilizado anteriormente. Entonces, las permutaciones son todas relativas a los conjuntos de datos de prueba originales que generamos.

A continuación, encontramos la correlación máxima:

cor_values = sapply(n_iterations, '[[', 'cor')
n_iterations[[which.max(cor_values)]]
$vec
  [1] 0.89799621 0.67532970 0.46456652 0.75948043 0.30868672 0.83244284
  [7] 0.86719899 0.55623242 0.63877904 0.73996913 0.71901990 0.85240338
 [13] 0.81109473 0.52571331 0.82931652 0.60854785 0.19701975 0.26986325
 [19] 0.58667068 0.52140366 0.40796732 0.22736077 0.74445997 0.40125309
 [25] 0.89193212 0.52031017 0.92285322 0.91790830 0.91780357 0.49734797
 [31] 0.07729027 0.11796154 0.69356590 0.95770015 0.74418811 0.43377516
 [37] 0.55981826 0.93903143 0.30047034 0.84776853 0.32203732 0.25235660
 [43] 0.79072381 0.58727094 0.99006038 0.01581969 0.41357545 0.52590383
 [49] 0.27980561 0.50573325 0.92423638 0.11261002 0.89840404 0.15345334
 [55] 0.61397396 0.27077191 0.12056045 0.45862163 0.18885955 0.77785348
 [61] 0.23440845 0.05244047 0.25931398 0.57217746 0.35554970 0.34622648
 [67] 0.21050067 0.08890968 0.84003379 0.95114398 0.83349287 0.82437619
 [73] 0.46724807 0.02641979 0.71740151 0.74439384 0.14830718 0.82685046
 [79] 0.33821824 0.71627101 0.77182339 0.72147500 0.08801319 0.08626022
 [85] 0.87332019 0.34303707 0.45393971 0.47871491 0.29254936 0.08939812
 [91] 0.35698405 0.67369873 0.27087693 0.78854984 0.87959539 0.22376395
 [97] 0.02674156 0.07077250 0.57461506 0.40616274

$cor
[1] 0.3166681

... o encuentre el valor más cercano a una correlación de 0.2:

n_iterations[[which.min(abs(cor_values - 0.2))]]
$vec
  [1] 0.02641979 0.49734797 0.32203732 0.95770015 0.82931652 0.52571331
  [7] 0.25931398 0.30047034 0.55981826 0.08801319 0.29254936 0.23440845
 [13] 0.12056045 0.89799621 0.57461506 0.99006038 0.27077191 0.08626022
 [19] 0.14830718 0.45393971 0.22376395 0.89840404 0.08890968 0.15345334
 [25] 0.87332019 0.92285322 0.50573325 0.40796732 0.91780357 0.57217746
 [31] 0.52590383 0.84003379 0.52031017 0.67532970 0.83244284 0.95114398
 [37] 0.81109473 0.35554970 0.92423638 0.83349287 0.34622648 0.18885955
 [43] 0.61397396 0.89193212 0.74445997 0.46724807 0.72147500 0.33821824
 [49] 0.71740151 0.75948043 0.52140366 0.69356590 0.41357545 0.21050067
 [55] 0.87959539 0.11796154 0.73996913 0.30868672 0.47871491 0.63877904
 [61] 0.22736077 0.40125309 0.02674156 0.26986325 0.43377516 0.07077250
 [67] 0.79072381 0.08939812 0.86719899 0.55623242 0.60854785 0.71627101
 [73] 0.40616274 0.35698405 0.67369873 0.82437619 0.27980561 0.77182339
 [79] 0.19701975 0.82685046 0.74418811 0.58667068 0.93903143 0.74439384
 [85] 0.46456652 0.85240338 0.34303707 0.45862163 0.91790830 0.84776853
 [91] 0.78854984 0.05244047 0.58727094 0.77785348 0.01581969 0.27087693
 [97] 0.07729027 0.71901990 0.25235660 0.11261002

$cor
[1] 0.2000199

Para obtener una mayor correlación, debe aumentar el número de iteraciones.

un problema más general: dada la variable ¿cómo generar las variables aleatorias con la matriz de correlación ?Y 2 , … , Y n$Y_1$$Y_2,\dots,Y_n$$R$

Solución:

1. obtener la descomposición cholesky de la matriz de correlación$CC^T=R$
2. crear vectores aleatorios independientes de la misma longitud queY $X_2,\dots,X_n$$Y_1$
3. Use como la primera columna y agregue los randoms generados$Y_1$
4. Y i Y $Y=CX$ , donde : los nuevos números aleatorios correlacionados según sea necesario, tenga en cuenta que no cambiará$Y_i$$Y_1$

Código de Python:

import numpy as np
import math
from scipy.linalg import toeplitz, cholesky
from statsmodels.stats.moment_helpers import cov2corr

# create the large correlation matrix R
p = 4
h = 2/p
v = np.linspace(1,-1+h,p)
R = cov2corr(toeplitz(v))

# create the first variable
T = 1000;
y = np.random.randn(T)

# generate p-1 correlated randoms
X = np.random.randn(T,p)
X[:,0] = y
C = cholesky(R)
Y = np.matmul(X,C)

# check that Y didn't change
print(np.max(np.abs(Y[:,0]-y)))

# check the correlation matrix
print(R)
print(np.corrcoef(np.transpose(Y)))

Prueba de salida:

0.0
[[ 1.   0.5  0.  -0.5]
 [ 0.5  1.   0.5  0. ]
 [ 0.   0.5  1.   0.5]
 [-0.5  0.   0.5  1. ]]
[[ 1.          0.50261766  0.02553882 -0.46259665]
 [ 0.50261766  1.          0.51162821  0.05748082]
 [ 0.02553882  0.51162821  1.          0.51403266]
 [-0.46259665  0.05748082  0.51403266  1.        ]]

Genere variables normales con la matriz de covarianza de MUESTREO como se indica

covsam <- function(nobs,covm, seed=1237) {; 
          library (expm);
          # nons=number of observations, covm = given covariance matrix ; 
          nvar <- ncol(covm); 
          tot <- nvar*nobs;
          dat <- matrix(rnorm(tot), ncol=nvar); 
          covmat <- cov(dat); 
          a2 <- sqrtm(solve(covmat)); 
          m2 <- sqrtm(covm);
          dat2 <- dat %*% a2 %*% m2 ; 
          rc <- cov(dat2);};
          cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3);
          cm; 
          res <- covsam(10,cm)  ;
          res;

Genere variables normales con la matriz de covarianza de POBLACIÓN como se indica

covpop <- function(nobs,covm, seed=1237) {; 
          library (expm); 
          # nons=number of observations, covm = given covariance matrix;
          nvar <- ncol(covm); 
          tot <- nvar*nobs;  
          dat <- matrix(rnorm(tot), ncol=nvar); 
          m2 <- sqrtm(covm);
          dat2 <- dat %*% m2;  
          rc <- cov(dat2); }; 
          cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3);
          cm; 
          res <- covpop(10,cm); 
          res

Simplemente cree un vector aleatorio y ordene hasta obtener el r deseado.