Uso de distribución uniforme para generar muestras aleatorias correlacionadas en R

[En preguntas recientes, estaba buscando generar vectores aleatorios en R , y quería compartir esa "investigación" como un Q&A independiente sobre un punto específico.]

La generación de datos aleatorios con correlación se puede hacer usando la descomposición de Cholesky de la matriz de correlación aquí , como se refleja en publicaciones anteriores aquí y aquí .$C = LL^{T}$

La pregunta que quiero abordar es cómo utilizar la distribución uniforme para generar números aleatorios correlacionados de diferentes distribuciones marginales en R .

Ya que la pregunta es

"cómo usar la distribución Uniforme para generar números aleatorios correlacionados de diferentes distribuciones marginales en "$\mathbb{R}$

y no solo las variaciones aleatorias normales, la respuesta anterior no produce simulaciones con la correlación prevista para un par arbitrario de distribuciones marginales en .$\mathbb{R}$

La razón es que, para la mayoría de los cdfs y , cuando donde denota el cdf normal estándar.$G_X$$G_Y$

A saber, aquí hay un contraejemplo con un Exp (1) y un Gamma (.2,1) como mi par de distribuciones marginales en .$\mathbb{R}$

library(mvtnorm)
#correlated normals with correlation 0.7
x=rmvnorm(1e4,mean=c(0,0),sigma=matrix(c(1,.7,.7,1),ncol=2),meth="chol")
cor(x[,1],x[,2])
  [1] 0.704503
y=pnorm(x) #correlated uniforms
cor(y[,1],y[,2])
  [1] 0.6860069
#correlated Exp(1) and Ga(.2,1)
cor(-log(1-y[,1]),qgamma(y[,2],shape=.2))
  [1] 0.5840085

Otro contraejemplo obvio es cuando es el cdf de Cauchy, en cuyo caso la correlación no está definida.$G_X$

Para dar una imagen más amplia, aquí hay un código R donde y son arbitrarios:$G_X$$G_Y$

etacor=function(rho=0,nsim=1e4,fx=qnorm,fy=qnorm){
  #generate a bivariate correlated normal sample
  x1=rnorm(nsim);x2=rnorm(nsim)
  if (length(rho)==1){
    y=pnorm(cbind(x1,rho*x1+sqrt((1-rho^2))*x2))
    return(cor(fx(y[,1]),fy(y[,2])))
    }
  coeur=rho
  rho2=sqrt(1-rho^2)
  for (t in 1:length(rho)){
     y=pnorm(cbind(x1,rho[t]*x1+rho2[t]*x2))
     coeur[t]=cor(fx(y[,1]),fy(y[,2]))}
  return(coeur)
  }

Jugar con diferentes cdfs me llevó a destacar este caso especial de una para y una distribución log-Normal para :$\chi^2_3$$G_X$$G_Y$

rhos=seq(-1,1,by=.01)
trancor=etacor(rho=rhos,fx=function(x){qchisq(x,df=3)},fy=qlnorm)
plot(rhos,trancor,ty="l",ylim=c(-1,1))
abline(a=0,b=1,lty=2)

que muestra qué tan lejos de la diagonal puede estar la correlación.

Una advertencia final Dadas dos distribuciones arbitrarias y , el rango de valores posibles de no es necesariamente . El problema puede no tener solución.$G_X$$G_Y$$\text{cor}(X,Y)$$(-1,1)$

Escribí el correlatepaquete. La gente dijo que es prometedor (digno de una publicación en el Journal of Statistical Software), pero nunca escribí el artículo porque elegí no seguir una carrera académica.

Creo que el correlatepaquete no mantenido todavía está en CRAN.

Cuando lo instales, puedes hacer lo siguiente:

require('correlate')
a <- rnorm(100)
b <- runif(100)
newdata <- correlate(cbind(a,b),0.5)

El resultado es que los nuevos datos tendrán una correlación de 0.5, sin cambiar las distribuciones univariadas de ay b(los mismos valores están allí, simplemente se mueven hasta que se alcanza la correlación multivariada 0.5.

Contestaré preguntas aquí, perdón por la falta de documentación.

1. Genere dos muestras de datos correlacionados a partir de una distribución aleatoria normal estándar siguiendo una correlación predeterminada .

   Como ejemplo, escojamos una correlación r = 0.7 y codifiquemos una matriz de correlación como:

   (C <- matrix(c(1,0.7,0.7,1), nrow = 2)) [,1] [,2] [1,] 1.0 0.7 [2,] 0.7 1.0

   Podemos usar mvtnormpara generar ahora estas dos muestras como un vector aleatorio bivariado:

   set.seed(0)

   SN <- rmvnorm(mean = c(0,0), sig = C, n = 1e5)resultando en dos componentes vectoriales distribuidos como ~ y con a . Ambos componentes se pueden extraer de la siguiente manera:$N(0, 1)$cor(SN[,1],SN[,2])= 0.6996197 ~ 0.7

   X1 <- SN[,1]; X2 <- SN[,2]

   Aquí está la trama con la línea de regresión superpuesta:

   

2. Use la Transformación integral de probabilidad aquí para obtener un vector aleatorio bivariado con distribuciones marginales ~$U(0, 1)$ y la misma correlación :

   U <- pnorm(SN)- entonces estamos alimentando pnormel SNvector para encontrar$erf(SN)$ (o $\Phi(SN)$) En el proceso, preservamos el cor(U[,1], U[,2]) = 0.6816123 ~ 0.7.

   Nuevamente podemos descomponer el vector U1 <- U[,1]; U2 <- U[,2]y producir un diagrama de dispersión con distribuciones marginales en los bordes, mostrando claramente su naturaleza uniforme:

   

3. Aplique el método de muestreo de transformación inversa aquí para obtener finalmente el bivector de puntos igualmente correlacionados que pertenezcan a cualquier familia de distribución que nos propongamos reproducir.

   A partir de aquí, podemos generar dos vectores distribuidos normalmente y con variaciones iguales o diferentes . Por ejemplo: Y1 <- qnorm(U1, mean = 8,sd = 10)y Y2 <- qnorm(U2, mean = -5, sd = 4), que mantendrá la correlación deseada, cor(Y1,Y2) = 0.6996197 ~ 0.7.

   O optar por diferentes distribuciones. Si las distribuciones elegidas son muy diferentes, la correlación puede no ser tan precisa. Por ejemplo, sigamos U1a$t$distribución con 3 df, y U2una exponencial con un$\lambda$= 1: Z1 <- qt(U1, df = 3)y Z2 <- qexp(U2, rate = 1)el cor(Z1,Z2) [1] 0.5941299 < 0.7. Aquí están los histogramas respectivos:

   

Aquí hay un ejemplo de código para todo el proceso y los márgenes normales:

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
}

A modo de comparación, he reunido una función basada en la descomposición de Cholesky:

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
}

Probar ambos métodos para generar correlacionados (por ejemplo, $r=0.7$) muestras distribuidas ~ $N(97,23)$ y $N(32,8)$obtenemos, configurando set.seed(99):

Usando el uniforme:

cor_samples(0.7, 1000, 97, 32, 23, 8)
           c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1, Y2))
mean Y1                                            96.5298821
mean Y2                                            32.1548306
SD Y1                                              22.8669448
SD Y2                                               8.1150780
cor(Y1,Y2)                                          0.7061308

y usando el Cholesky:

Cholesky_samples(0.7, 1000, 97, 32, 23, 8)
             c(mean(N_1), mean(N_2), sd(N_1), sd(N_2), cor(N_1, N_2))
mean N_1                                                   96.4457504
mean N_2                                                   31.9979675
SD N_1                                                     23.5255419
SD N_2                                                      8.1459100
cor(N_1,N_2)                                                0.7282176