Utilisation d'une distribution uniforme pour générer des échantillons aléatoires corrélés dans R

[Sur des questions récentes, je cherchais à générer des vecteurs aléatoires dans R , et je voulais partager cette "recherche" en tant que Q&A indépendante sur un point spécifique.]

La génération de données aléatoires avec corrélation peut être effectuée en utilisant la décomposition de Cholesky de la matrice de corrélation ici , comme reflété dans les articles précédents ici et ici .$C = LL^{T}$

La question que je veux aborder est comment utiliser la distribution uniforme pour générer des nombres aléatoires décorrélés de différentes distributions marginales en R .

La question étant

"comment utiliser la distribution uniforme pour générer des nombres aléatoires corrélés à partir de différentes distributions marginales dans "$\mathbb{R}$

et pas seulement des variables aléatoires normales, la réponse ci-dessus ne produit pas de simulations avec la corrélation voulue pour une paire arbitraire de distributions marginales dans .$\mathbb{R}$

La raison en est que, pour la plupart des cdfs et , lorsque où désigne le cdf normal standard.$G_X$$G_Y$

À savoir, voici un contre-exemple avec un Exp (1) et un Gamma (.2,1) comme ma paire de distributions marginales dans .$\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

Un autre contre-exemple évident est quand est le cdf de Cauchy, auquel cas la corrélation n'est pas définie.$G_X$

Pour donner une image plus large, voici un code R où et sont arbitraires:$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)
  }

Jouer avec différents cdfs m'a amené à distinguer ce cas particulier d'une pour et d'une distribution log-Normal pour :$\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)

ce qui montre à quelle distance de la diagonale la corrélation peut être.

Un dernier avertissement Étant donné deux distributions arbitraires et , la plage de valeurs possibles de n'est pas nécessairement . Le problème peut donc n'avoir aucune solution.$G_X$$G_Y$$\text{cor}(X,Y)$$(-1,1)$

J'ai écrit le correlatepaquet. Les gens ont dit que c'était prometteur (digne d'une publication dans le Journal of Statistical Software), mais je n'ai jamais écrit l'article pour cela parce que j'ai choisi de ne pas poursuivre une carrière universitaire.

Je crois que le correlatepackage non maintenu est toujours sur CRAN.

Lorsque vous l'installez, vous pouvez effectuer les opérations suivantes:

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

Le résultat est que les nouvelles données auront une corrélation de 0,5, sans changer les distributions univariées de aet b(les mêmes valeurs sont là, elles sont simplement déplacées jusqu'à ce que la corrélation multivariée 0,5 soit atteinte.

Je répondrai aux questions ici, désolé pour le manque de documentation.

1. Générez deux échantillons de données corrélées à partir d'une distribution aléatoire normale standard suite à une corrélation prédéterminée .

   Par exemple, prenons une corrélation r = 0,7 et codons une matrice de corrélation telle que:

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

   Nous pouvons utiliser mvtnormpour générer maintenant ces deux échantillons comme un vecteur aléatoire bivarié:

   set.seed(0)

   SN <- rmvnorm(mean = c(0,0), sig = C, n = 1e5)résultant en deux composantes vectorielles distribuées comme ~ et avec a . Les deux composants peuvent être extraits comme suit:$N(0, 1)$cor(SN[,1],SN[,2])= 0.6996197 ~ 0.7

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

   Voici l'intrigue avec la ligne de régression qui se chevauche:

   

2. Utilisez la transformation intégrale de probabilité ici pour obtenir un vecteur aléatoire bivarié avec des distributions marginales ~ et la même corrélation :$U(0, 1)$

   U <- pnorm(SN)- donc nous alimentons pnormle SNvecteur pour trouver (ou ). Ce faisant, nous préservons le .$erf(SN)$$\Phi(SN)$cor(U[,1], U[,2]) = 0.6816123 ~ 0.7

   Encore une fois, nous pouvons décomposer le vecteur U1 <- U[,1]; U2 <- U[,2]et produire un diagramme de dispersion avec des distributions marginales sur les bords, montrant clairement leur nature uniforme:

   

3. Appliquez ici la méthode d'échantillonnage par transformée inverse pour finalement obtenir le bivecteur de points également corrélés appartenant à la famille de distribution que nous nous proposons de reproduire.

   De là, nous pouvons simplement générer deux vecteurs distribués normalement et avec des variances égales ou différentes . Par exemple: Y1 <- qnorm(U1, mean = 8,sd = 10)et Y2 <- qnorm(U2, mean = -5, sd = 4)qui maintiendra la corrélation désirée, cor(Y1,Y2) = 0.6996197 ~ 0.7.

   Ou optez pour différentes distributions. Si les distributions choisies sont très différentes, la corrélation peut ne pas être aussi précise. Par exemple, U1suivons une distribution avec 3 df, et une exponentielle avec a = 1: et The . Voici les histogrammes respectifs:$t$U2$\lambda$Z1 <- qt(U1, df = 3)Z2 <- qexp(U2, rate = 1)cor(Z1,Z2) [1] 0.5941299 < 0.7

   

Voici un exemple de code pour tout le processus et les marginaux normaux:

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
}

À titre de comparaison, j'ai mis en place une fonction basée sur la décomposition 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
}

En essayant les deux méthodes pour générer des échantillons corrélés (disons ) distribués ~ et nous obtenons, en définissant :$r=0.7$$N(97,23)$$N(32,8)$set.seed(99)

Utilisation de l'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

et utilisation du 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