2012-05-10 11 views
10

複数の二項分布から相関した乱数を生成する方法を見つけようとしています。二項分布から相関した乱数を生成するR

私は通常のディストリビューション(mvrnormを使って)でそれを行う方法を知っていますが、私はバイナリのものに適用できる関数を見つけませんでした。

+0

パッケージ 'bindata'を使うことができます。ここでうまくデモします:https://stat.ethz.ch/pipermail/r-help/2007-July/135575.html。 (そのリンクは、 'Rの相関二項変数をシミュレートするためのGoogleの検索で返された最初のページにありました...) –

+0

ありがとうございました。ジョシュですが、バイナリデータではなく2項データが必要です。 – Arnaud

+1

@Arnaud - 私は今朝カフェインや覚せい剤を一切持っていないが、唯一の許容値が「はい/いいえ」、「合格/不合格」、「真/ FALSE "、つまりバイナリ?それは[ウィキペディアも考えているようだ](http://en.wikipedia.org/wiki/Binomial_distribution) – Chase

答えて

11

の場合はそうではありません。ここでは1つの簡単な例です:

library(copula) 

tmp <- normalCopula(0.75, dim=2) 
x <- rcopula(tmp, 1000) 
x2 <- cbind(qbinom(x[,1], 10, 0.5), qbinom(x[,2], 15, 0.7)) 

x2が相関している2二項変数を表す2列の行列です。

+0

うわー、素敵で甘い! – Arnaud

+0

もう一度Greg ... Rの助けを借りてあなたの助けを借りて、もう一度私を救ってください! – Arnaud

+4

これは興味深い考えですが、目的の相関を持つ変数を返しません。 (例えば、私は上記のコード100回分のサンプル相関係数を計算した:平均相関は0.724であり、相関係数のうちのわずか5つは0.75より大きい)。 –

9

各試行でn回の試行と成功確率pを持つ2項変数は、それぞれ成功確率pが であるn回のベルヌーイ試行の合計と見ることができます。

同様に、所望の相関rを有するベルヌーイ変量の対を合計することによって、相関二項変量の対を で構成することができます。

require(bindata) 

# Parameters of joint distribution 
size <- 20 
p1 <- 0.5 
p2 <- 0.3 
rho<- 0.2 

# Create one pair of correlated binomial values 
trials <- rmvbin(size, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) 
colSums(trials) 

# A function to create n correlated pairs 
rmvBinomial <- function(n, size, p1, p2, rho) { 
    X <- replicate(n, { 
      colSums(rmvbin(size, c(p1,p2), bincorr=(1-rho)*diag(2)+rho)) 
     }) 
    t(X) 
} 
# Try it out, creating 1000 pairs 
X <- rmvBinomial(1000, size=size, p1=p1, p2=p2, rho=rho) 
#  cor(X[,1], X[,2]) 
# [1] 0.1935928 # (In ~8 trials, sample correlations ranged between 0.15 & 0.25) 

それは希望の相関係数を共有多くの異なる結合分布があることに注意することが重要です。 rmvBinomial()のシミュレーション方法では、その1つが生成されますが、それが適切かどうかは、データを生成するプロセスに依存します。

同様の問題(次に詳細にアイデアを説明 に進む)にthis R-help answerに述べたように、(与えられた平均及び分散)二変量正常を一意相関係数によって定義されている間

、これはあなたが二項変数にそれらを変換するためにqbinom機能を使用し、その後、copulaパッケージを使用して相関ユニフォームを生成することができ、二変量二項

+0

ありがとうジョシュ。より多くの数の二項分布を可能にするようにスクリプトを修正しました。しかし、http://stat.ethz.ch/pipermail/r-help/2007-July/135575.htmlに示されているように、rhoは限界の確率のある関数によって下と上に境界が定められている(関数はrho = 0.8で失敗する)。奇数の比率の使用は解決策だと思われますが... 2つ以上の分布に対して有効な2進数の相関にオッズ比を変換する機能を一般化する方法を知っていますか? – Arnaud

+0

@Josh私は関連する質問をしましたが、おそらくあなたはそれを見たいでしょうか? https://stackoverflow.com/questions/47006881/how-to-generate-a-binomial-vector-of-n-correlated-items – jaySf

関連する問題