r - それらのペア間の特定の相関関係を持つ確率変数の生成:

Question

2 つの連続確率変数 (量的特性、それぞれが正常) と 2 つのバイナリ確率変数 (バイナリ特性) を生成Q1しQ2、Z1それらZ2のすべての可能なペア間のペアワイズ相関を指定します。言う

(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でそのようなデータを生成するのを手伝ってください.

Answer 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

相関値をシミュレートする (最初の 2 つの定量的値、最後の 2 つをバイナリに変換)

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

以前の値よりもはるかに優れているようには見えません。考えられる 2 つの問題 (読者の演習として残しておきます) (1) 選択した周辺分布では実現不可能な一連の相関を指定した、または (2) 目的関数曲面の誤差が方法 - より良く行うには、より多くの複製を平均化する必要があります (つまり、 を増やしnます)。