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

71

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]

Felix S
fuente
2
Consulte esta pregunta relacionada stats.stackexchange.com/questions/13382/… que aborda directamente su pregunta (al menos el lado teórico).
Macro
La siguiente Q también está fuertemente relacionada y será de interés: Cómo generar números aleatorios correlacionados (dado significa varianzas y grado de correlación) .
gung - Restablecer Monica

Respuestas:

56

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 θxrθ

  1. obtener el vector fijo y un vector aleatoriox 2x1x2
  2. centrar ambos vectores (media 0), dando vectores , ˙ x 2x˙1x˙2
  3. hacer ortogonal a (proyección sobre el subespacio ortogonal), dando ˙ x 1 ˙ x 2x˙2x˙1x˙2
  4. escala y a la longitud 1, dando y ˙ x 2 ˉ x 1 ˉ x 2x˙1x˙2x¯1x¯2
  5. ˉ x 1θ ˉ x 1rx1x¯2+(1/tan(θ))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.x¯1θx¯1rx1

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

ingrese la descripción de la imagen aquí

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

lince
fuente
Estaba tratando de reescribir el código en la sintaxis de SPSS. Me tropiezo con tu descomposición QR que devuelve una columna de 20x1. En SPSS tengo la ortonormalización de Gram-Schmidt (que también es una descomposición QR) pero no puedo replicar la columna Q resultante. ¿Puedes masticarme tu acción QR por favor? O indique alguna solución para obtener la proyección. Gracias.
ttnphns
@caracal, P <- X %*% solve(t(X) %*% X) %*% t(X)no produce r = 0.6, así que esa no es la solución . Todavía estoy confundido. (Estaré encantado de imitar su expresión Q <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))en SPSS pero no sé cómo.)
ttnphns
@ttnphns Perdón por la confusión, mi comentario fue para el caso general. Aplicandolo a la situación en el ejemplo: Obtener la matriz de proyección mediante descomposición QR es solo para estabilidad numérica. Puede obtener la matriz de proyección como si las columnas de la matriz abarcan el subespacio . En R, puede escribir aquí porque el subespacio está abarcado por la primera columna de . La matriz para la proyección sobre el complemento ortogonal es entonces IP. XP=X(XX)1XXXctr[ , 1] %*% solve(t(Xctr[ , 1]) %*% Xctr[ , 1]) %*% t(Xctr[ , 1])Xctr
caracal
44
¿Alguien podría aclarar cómo realizar algo similar para más de solo dos muestras? Digamos, si quisiera 3 muestras correlacionadas por pares por rho, ¿cómo puedo transformar esta solución para lograr eso?
Andre Terra
para el caso límite, rho=1me pareció útil hacer algo como esto: de lo if (isTRUE(all.equal(rho, 1))) rho <- 1-10*.Machine$double.epscontrario, estaba recibiendo NaNs
PatrickT el
19

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 ρ YXYYXYYXYYρY

XY;ρ=ρSD(Y)Y+1ρ2SD(Y)Y.

(" " representa cualquier cálculo proporcional a una desviación estándar).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 ) YY50XY;ρYX=(1,2,,50)Y

Figura

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 YYXY1,Y2,,Yk;ρ1,ρ2,,ρkYiYiXYiYY

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

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))
whuber
fuente
De hecho, esta es una buena solución. Sin embargo, no pude expandirlo yo mismo a múltiples variables (las variables fijas, en su respuesta). , tu aseguras. ¿Puedes demostrarlo? Por favor, ¿con un código anotado legible por un usuario que no sea R? YBTW, this method readily generalizes to more... Just use ordinary least squares... and form a suitable linear combination
ttnphns
1
@ttnphns lo he hecho.
whuber
1
¡Muchas gracias! Ya veo, y he codificado su enfoque hoy en SPSS para mí. Muy buena propuesta tuya. Nunca pensé en la noción de base dual como aplicable para resolver la tarea.
ttnphns
¿Es posible utilizar un enfoque similar para obtener un vector distribuido uniformemente? Es decir, tengo un vector existente xy quiero generar un nuevo vector ycorrelacionado xpero también quiero que el yvector se distribuya uniformemente.
Skumin
@Skumin Considere usar una cópula para eso para que pueda controlar la relación entre los dos vectores.
whuber
6

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 .ρρρ

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:ρ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).

Felix S
fuente
Parece que esto es solo una solución aproximada, es decir, la correlación empírica no es exactamente igual a . ¿O me estoy perdiendo algo? ρ
caracal
1
Es fácil demostrar que, excepto por esa "pequeña corrección a rho" (cuyo propósito en este contexto me elude), esto es exactamente lo mismo que lo sugerido anteriormente por ttnphns. El método se basa simplemente en la descomposición de Choleski de la matriz de correlación para obtener la matriz de transformación deseada. Ver, por ejemplo: en.wikipedia.org/wiki/… . Y sí, esto solo le dará dos vectores cuya correlación de población es igual a rho.
Wolfgang
La "pequeña corrección a rho" estaba en la publicación original y se describe aquí . En realidad, realmente no lo entiendo; pero una investigación de 50000 correlaciones simuladas con rho = .3 muestra que sin la "pequeña corrección" se produce un promedio de r's de .299, mientras que con la corrección se obtiene un promedio de .312 (que es el valor del rho corregido) producido. Por lo tanto, eliminé esa parte de la función.
Felix S
Sé que esto es antiguo, pero también quiero señalar que este método no funcionará para matrices de correlación definidas no positivas. Por ejemplo, una correlación de -1.
zzk
1
Gracias; Me di cuenta de que si x1 no ha sido estandarizada media = 0, SD = 1, y prefiere no cambiar la escala de él, tendrá que modificar la línea: X2 <- mar.fun(n)a X2 <- mar.fun(n,mean(x),sd(x))obtener la correlación deseada entre x1 y x2
Dave M
6

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 (siXYXrXrY=rX+EE0sd=1r2XYrXYX es de lo normal) con .ρ=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.rEXEXYX1,X2,X3,...

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.XrYYrY


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.X1,X2,...,XmYYr1,r2,...,rmX

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).YXYY

  1. Convierta el objetivo s en sumas de productos cruzados multiplicándolos por : . ( es un índice variable ).rdf=n1Sj=rjdfjX

  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 = .dfYXdf

  3. Calcular los coeficientes de predicción de regressional por s de acuerdo con el objetivo s: .YXrb=(XX)1S

  4. Calcule los valores pronosticados para : .YY^=Xb

  5. Calcule los residuos .E=YY^

  6. Calcule la suma de cuadrados (objetivo) necesaria para los residuos: .SSS=dfSSY^

  7. (Comience a repetir.) Calcule las sumas observadas de productos cruzados entre actual y cada :EXjCj=i=1nEiXij

  8. Corrija los valores de con el objetivo de acercar todos los s a ( es un índice de caso):EC0i

    Ei[corrected]=Eij=1mCjXijnj=1mXij2

    (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:E0 EC

    Ei[corrected]=Eij=1mCjXij3i=1nXij2j=1mXij2

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

  9. Traiga a su valor objetivo:SSEEi[corrected]=EiSSS/SSE

    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.)mrSSSn

  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 + ECErYY[corrected]=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).rYr

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


Y X1 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 ).YX

ttnphns
fuente
1
Gracias por tu respuesta. Esa es una solución empírica / iterativa en la que también estaba pensando. Sin embargo, para mis simulaciones, necesito una solución más analítica sin un procedimiento de ajuste costoso. Afortunadamente, acabo de encontrar una solución que publicaré en breve ...
Felix S
Esto funciona para generar normales bivariadas, pero no funciona para una distribución arbitraria (o cualquier distribución no 'aditiva')
Macro
1
No veo por qué propone una iteración cuando puede producir todo el cono de soluciones directamente. ¿Hay algún propósito especial para este enfoque?
whuber
1
Re su última edición: dado que proporciono una fórmula simple para todas las soluciones, se puede lograr cualquier objetivo deseado, como "mayor homocedasticidad", minimizando una función objetivo adecuada sobre el conjunto de todas las soluciones. El enfoque es totalmente general. Al extender la variable (o variables) a una base ortogonal y explotar la invariancia de escala de la correlación, el problema se convierte en optimizar una función definida en una esfera en un espacio euclidiano. Y
whuber
1
@whuber, tu comentario es lo que estaba esperando; En realidad, mi respuesta (sobre la heterocedasticidad, a la que me relaciono) fue un desafío para usted: tal vez es una invitación a publicar su solución, tan minuciosa y brillante como lo hace habitualmente.
ttnphns
4

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.

Paul Hiemstra
fuente
2

un problema más general: dada la variable ¿cómo generar las variables aleatorias con la matriz de correlación ?Y 2 , , Y n RY1Y2,,YnR

Solución:

  1. obtener la descomposición cholesky de la matriz de correlaciónCCT=R
  2. crear vectores aleatorios independientes de la misma longitud queY 1X2,,XnY1
  3. Use como la primera columna y agregue los randoms generadosY1
  4. Y i Y 1Y=CX , donde : los nuevos números aleatorios correlacionados según sea necesario, tenga en cuenta que no cambiaráYiY1

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.        ]]
Aksakal
fuente
¿Podría aclarar qué significa "no que no cambiará"? Y1
whuber
@whuber fue un error tipográfico
Aksakal
0

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
usuario3635627
fuente
2
¡Necesita aprender a formatear el código en la respuesta! Hay una opción específica para marcar texto como fragmentos de código, ¡úselo!
kjetil b halvorsen
-6

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

Adán
fuente
¿En qué situaciones sería preferible esto a las soluciones anteriores?
Andy W
Una situación en la que un usuario quiere una respuesta simple. Leí una pregunta similar en el foro r, y es la respuesta que se dio.
Adam
3
Desafortunadamente, esta solución no solo es computacionalmente ineficiente y aproximada, a menudo fallará por completo a menos que se aplique primero algún análisis para determinar una distribución apropiada para el "vector aleatorio". Creo que tiene mérito la idea subyacente de simplemente arrojar algunos números aleatorios al problema y permutarlos aleatoriamenteno "ordenarlos") hasta que se obtenga una aproximada (porque esto es rápido y fácil de programar), pero esa idea no se expresa claramente en esta breve respuesta. r
whuber
3
Si esta respuesta se dio en el foro de r-help, sospecho que fue (a) irónico (es decir, como una broma) u (b) ofrecido por alguien que no es estadísticamente muy sofisticado. Para decirlo de manera más sucinta, esta es una respuesta pobre a la pregunta. -1
gung - Restablece a Monica