【问题标题】:Generating Random Variables with given correlations between pairs of them:生成具有给定相关性的随机变量:
【发布时间】:2014-07-13 19:53:46
【问题描述】:

我想生成 2 个连续随机变量 Q1Q2(定量特征,每个都是正常的)和 2 个二元随机变量 Z1Z2(二元特征),在所有可能对之间具有给定的成对相关性其中。 说

(Q1,Q2):0.23 
(Q1,Z1):0.55 
(Q1,Z2):0.45 
(Q2,Z1):0.4 
(Q2,Z2):0.5 
(Z1,Z2):0.47 

请帮助我在 R 中生成此类数据。

【问题讨论】:

  • 为什么会涨票? 1)缺乏关键信息; 2) 含糊不清; 3) 没有表现出努力; 4)“请为我做点事;” 5) 叹息
  • 对我来说 (1) 和 (2) 似乎不是真的,但 (3) 和 (4) 是。我不知道有一种简单/简单的方法可以做到这一点,但它看起来确实很好。大概是因为人们对答案感兴趣。
  • 您可以轻松地独立于定性特征的平均值指定相关性(例如,使它们成为多元正态),但您可能不得不说一下二元随机变量的平均概率。跨度>
  • 如果 op 不想要 mvn 怎么办?你不认为这是关键信息吗?
  • 实际上,我无法表现出任何努力,因为我不知道该怎么做。我可以生成相关的正态变量和相关的伯努利变量,但这一个超出了我的能力范围。另外,如果有人帮助我编写代码或建议可能的某种方式,我将不胜感激。

标签: r statistics simulation genetics


【解决方案1】:

这很粗略,但可能会让你朝着正确的方向开始。

library(copula)

options(digits=3)
probs <- c(0.5,0.5)
corrs <- c(0.23,0.55,0.45,0.4,0.5,0.47)  ## lower triangle

模拟相关值(前两个定量,后两个转换为二进制)

sim <- function(n,probs,corrs) {
    tmp <- normalCopula( corrs, dim=4 , "un")
    getSigma(tmp) ## test
    x <- rCopula(1000, tmp)
    x2 <- x
    x2[,3:4] <- qbinom(x[,3:4],size=1,prob=rep(probs,each=nrow(x)))
    x2
}

测试观察到的和目标相关性之间的 SSQ 距离:

objfun <- function(corrs,targetcorrs,probs,n=1000) {
    cc <- try(cor(sim(n,probs,corrs)),silent=TRUE)
    if (is(cc,"try-error")) return(NA)
    sum((cc[lower.tri(cc)]-targetcorrs)^2)
}

当输入 corrs=target 时,看看事情有多糟糕:

cc0 <- cor(sim(1000,probs=probs,corrs=corrs))
cc0[lower.tri(cc0)]
corrs
objfun(corrs,corrs,probs=probs) ## 0.112

现在尝试优化。

opt1 <- optim(fn=objfun,
              par=corrs,
              targetcorrs=corrs,probs=c(0.5,0.5))
opt1$value     ## 0.0208

在 501 次迭代后停止,“超出最大迭代次数”。这永远不会很好地工作,因为我们正在尝试在随机目标函数上使用确定性爬山算法......

cc1 <- cor(sim(1000,probs=c(0.5,0.5),corrs=opt1$par))
cc1[lower.tri(cc1)]
corrs

也许试试模拟退火?

opt2 <- optim(fn=objfun,
              par=corrs,
              targetcorrs=corrs,probs=c(0.5,0.5),
              method="SANN")

它似乎并没有比之前的值好多少。两个可能的问题(留给读者作为练习)(1)我们指定了一组与我们选择的边际分布不可行的相关性,或者(2)目标函数表面的误差在方式——为了做得更好,我们必须对更多的重复进行平均(即增加n)。

【讨论】:

    猜你喜欢
    • 2015-02-27
    • 1970-01-01
    • 2020-02-19
    • 1970-01-01
    • 1970-01-01
    • 2020-01-30
    • 1970-01-01
    • 1970-01-01
    • 2015-07-20
    相关资源
    最近更新 更多