在 R 中使用均匀分布生成相关随机样本

机器算法验证 r 相关性 采样 随机变量 随机生成
2022-03-14 06:03:18

[关于最近的问题,我正在研究在R中生成随机向量,我想将这项“研究”作为关于特定点的独立问答分享。]

可以在这里使用相关矩阵的 Cholesky 分解来生成具有相关性的随机数据,这反映在此处此处的先前帖子中。C=LLT

我要解决的问题是如何使用均匀分布从R中的不同边缘分布生成相关随机数。

3个回答

既然问题是

“如何使用均匀分布从中的不同边际分布生成相关随机数”R

不仅是正态随机变量,上述答案不会产生与中任意一对边际分布的预期相关性的模拟R

原因是,对于大多数 cdfs其中表示标准法线 cdf。GXGY

cor(X,Y)cor(GX1(Φ(X),GY1(Φ(Y)),
(X,Y)N2(0,Σ),
Φ

也就是说,这是一个反例,其中 Exp(1) 和 Gamma(.2,1) 作为我在中的一对边际分布。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

另一个明显的反例是当是 Cauchy cdf 时,在这种情况下没有定义相关性。GX

为了给出更广泛的图景,这是一个 R 代码,其中都是任意的:GXGY

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

在此处输入图像描述

玩弄不同的 cdfs 让我挑出了 G_X 的分布和的对数正态分布的这种特殊情况χ32GXGY

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)

这表明相关性可以离对角线多远。

最后警告给定两个任意分布的可能范围不一定是 因此问题可能没有解决方案。GXGYcor(X,Y)(1,1)

我写了correlate包裹。人们说它很有前途(值得在 Journal of Statistical Software 上发表),但我从来没有为它写过论文,因为我选择不从事学术生涯。

我相信未维护的correlate软件包仍在 CRAN 上。

安装时,您可以执行以下操作:

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

结果是 newdata 将具有 0.5 的相关性,而不会改变 和 的单变量分布ab存在相同的值,它们只是四处移动,直到达到多变量 0.5 相关性。

我会在这里回答问题,抱歉缺少文档。

  1. 根据预先确定的相关性,从标准正态随机分布中生成两个相关数据样本

    例如,让我们选择一个相关性r = 0.7,并编写一个相关性矩阵,例如:

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

    我们现在可以使用mvtnorm将这两个样本生成为二元随机向量:

    set.seed(0)

    SN <- rmvnorm(mean = c(0,0), sig = C, n = 1e5)导致两个向量分量分布为 ~N(0,1)并带有cor(SN[,1],SN[,2])= 0.6996197 ~ 0.7. 这两个组件都可以按如下方式提取:

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

    这是具有重叠回归线的图:

  2. 在这里使用概率积分变换 来获得具有边际分布的二元随机向量~U(0,1) 和相同的相关性

    U <- pnorm(SN)- 所以我们正在输入pnorm向量SN来查找erf(SN)(或者Φ(SN))。在此过程中,我们保留cor(U[,1], U[,2]) = 0.6816123 ~ 0.7.

    再次,我们可以分解向量U1 <- U[,1]; U2 <- U[,2]并生成边缘边缘分布的散点图,清楚地显示它们的均匀性质:

  3. 在这里应用逆变换采样方法,最终获得属于我们要重现的任何分布族的等相关点的双向量。

    从这里我们可以生成两个正态分布且方差相等或不同的向量。例如:Y1 <- qnorm(U1, mean = 8,sd = 10)Y2 <- qnorm(U2, mean = -5, sd = 4),这将保持所需的相关性,cor(Y1,Y2) = 0.6996197 ~ 0.7

    或者选择不同的分布。如果选择的分布非常不同,则相关性可能不那么精确。例如,让我们U1关注一个t具有 3 df 的分布和U2具有 a 的指数分布λ=1 :Z1 <- qt(U1, df = 3)Z2 <- qexp(U2, rate = 1). cor(Z1,Z2) [1] 0.5941299 < 0.7以下是各自的直方图:

这是整个过程和正常边缘的代码示例:

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
}

为了比较,我整理了一个基于 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
}

尝试这两种方法来生成相关(例如,r=0.7) 样本分发~N(97,23)N(32,8)我们得到,设置set.seed(99)

使用制服:

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

并使用 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