Utilizando distribución uniforme para generar muestras aleatorias correlacionadas r

Question

[En los últimos preguntas que yo estaba buscando en la generación aleatoria de los vectores en R, y quería compartir esa "investigación" como un independiente de Q&A sobre un punto específico.]

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

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

Answer 1

Ya que la pregunta es

"cómo utilizar la distribución Uniforme para generar correlacionado al azar los números de diferentes distribuciones marginales en $\mathbb{R}$"

y no sólo aleatoria normal varia, a la respuesta anterior, no produce simulaciones con la intención de correlación para cualquier par de distribuciones marginales en $\mathbb{R}$.

La razón es que, para la mayoría de los cdfs $G_X$ y $G_Y$,$$\text{cor}(X,Y)\ne\text{cor}(G_X^{-1}(\Phi(X),G_Y^{-1}(\Phi(Y)),$$when$$(X,Y)\sim\text{N}_2(0,\Sigma),$$where $\Phi$ denotes the standard normal cdf.

To wit, here is a counter-example with an Exp(1) and a Gamma(.2,1) as my pair of marginal distributions in $\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

Another obvious counter-example is when $G_X$ is the Cauchy cdf, in which case the correlation is not defined.

To give a broader picture, here is an R code where both $G_X$ and $G_Y$ are arbitrary:

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)
  }

Playing around with different cdfs led me to single out this special case of a $\chi^2_3$ distribution for $G_X$ and a log-Normal distribution for $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)

la que se muestra cómo la medida de la diagonal de la correlación puede ser.

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

Answer 2

Escribí el correlate paquete. La gente dice que es prometedor (digno de una publicación en el Journal of Statistical Software), pero nunca me escribió el papel para él, porque he decidido no seguir adelante con una carrera académica.

Yo creo que el no se mantienen correlate paquete está todavía en CRAN.

Cuando lo instale, puede hacer lo siguiente:

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

El resultado es que nuevosdatos tendrá una correlación de 0.5, sin cambiar la univariado de las distribuciones de a y b (mismos valores están ahí, que acaban de llegar se trasladó hasta el multivariante 0.5 correlación ha sido alcanzado.

Te voy a responder en preguntas aquí, lo siento por la falta de documentación.

Answer 3

1. Generar dos muestras de datos correlacionados de una normal estándar distribución aleatoria siguiendo una determinada correlación.

   Como ejemplo, vamos a elegir una correlación de r = 0.7, y el código de una matriz de correlación, tales 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 mvtnorm a generar ahora estas dos muestras como un bivariante vector aleatorio:

   set.seed(0)

   SN <- rmvnorm(mean = c(0,0), sig = C, n = 1e5) resultante de dos vectores componentes distribuidos como ~ $N(0, 1)$ y con un cor(SN[,1],SN[,2])= 0.6996197 ~ 0.7. Ambos componentes pueden ser llevadas de la siguiente manera:

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

   Aquí está la parcela con la superposición de la línea de regresión:

   

2. El uso de la Probabilidad de Transformación Integral aquí para obtener una bivariante vector aleatorio con distribuciones marginales ~ $U(0, 1)$ e la misma correlación:

   U <- pnorm(SN) - por lo que estamos en la alimentación de pnorm el SN vector para encontrar $erf(SN)$ (o $\Phi(SN)$). En el proceso, podemos preservar la cor(U[,1], U[,2]) = 0.6816123 ~ 0.7 .

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

   

3. Aplicar la inversa de la transformación método de muestreo aquí para obtener finalmente la bivector de igual correlaciona los puntos pertenecientes a cualquiera de distribución de la familia nos dispusimos a reproducir.

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

   O bien optar por diferentes distribuciones. Si las distribuciones elegidas son muy diferentes, la correlación puede no ser tan precisa. Por ejemplo, veamos U1 seguir $t$ distribución con 3 d.f., y U2 una 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 respectivos histogramas:

   

Aquí es un ejemplo de código para el proceso completo y normal marginales:

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
}

Para la comparación, he reunido una función que se basa 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
}

Tratando de ambos métodos para generar correlacionados (por ejemplo, $r=0.7$) muestras distribuidas ~ $N(97,23)$ $N(32,8)$ somos, estableciendo 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 el Uso de la 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