2010-02-17 40 views
5

假设我有一个二元离散分布,即对于i = 1,...,n和j = 1的概率值表P(X = i,Y = j)。 ..m。如何从这样的分布生成一个随机样本(X_k,Y_k),k = 1,... N?也许有一个现成的R功能,如:给定二元离散分布的随机样本

sample(100,prob=biprob) 

其中biprob是2维矩阵?

一个直观的方法是采样如下。假设我们有一个data.frame

dt=data.frame(X=x,Y=y,P=pij) 

其中X和Y来自

expand.grid(x=1:n,y=1:m) 

和伊斯兰圣战是P(X = I,Y = j)的。

然后我们得到我们的样本大小为N,通过以下方式(XS,YS):

set.seed(1000) 
Xs <- sample(dt$X,size=N,prob=dt$P) 
set.seed(1000) 
Ys <- sample(dt$Y,size=N,prob=dt$P) 

我用set.seed()来模拟 “bivariateness”。直觉上我应该得到类似于我需要的东西。我不确定这是否正确。因此,这个问题:)

另一种方法是使用吉布斯抽样,边际分布很容易计算。

我尝试了谷歌搜索,但没有真正相关出现。

回答

7

你快到了。假设你有数据帧dt与x,y和pij值,只需对行进行采样!

dt <- expand.grid(X=1:3, Y=1:2) 
dt$p <- runif(6) 
dt$p <- dt$p/sum(dt$p) # get fake probabilities 
idx <- sample(1:nrow(dt), size=8, replace=TRUE, prob=dt$p) 
sampled.x <- dt$X[idx] 
sampled.y <- dt$Y[idx] 
+0

再仔细读这篇文章,这是相同的解决方案,我建议。采样行可能比组合rmultinom和哪个更清晰。关键是要认识到行和列只是符号。 – Tristan

+0

是符号是关键。二元离散分布与单变量离散分布相同,符号改变。我选择Anika的答案是正确的,但仅仅因为代码更简单:) Tristan给出了更好的理论解释。 – mpiktas

+0

+1为好例子 – andi

7

我不清楚为什么你应该关心它是二元的。概率总和为1,结果是离散的,所以你只是从categorical distribution抽样。唯一的区别是您使用行和列而不是单个位置对观察值进行索引。这只是表示法。

在R中,您可以通过重新定型数据和从分类分布中抽样来轻松地从分销中抽样。可以使用rmultinom并使用which来选择索引,或者如Aniko所示,使用sample对重新整形数据的行进行采样,从而可以从分类中进行采样。一些簿记可以照顾你的确切情况。

这里有一个解决方案:

library(reshape) 

# Reshape data to long format. 
data <- matrix(data = c(.25,.5,.1,.4), nrow=2, ncol=2) 
pmatrix <- melt(data) 

# Sample categorical n times. 
rcat <- function(n, pmatrix) { 
    rows <- which(rmultinom(n,1,pmatrix$value)==1, arr.ind=TRUE)[,'row'] 
    indices <- pmatrix[rows, c('X1','X2')] 
    colnames(indices) <- c('i','j') 
    rownames(indices) <- seq(1,nrow(indices)) 
    return(indices) 
} 

rcat(3,pmatrix) 

这将返回3个随机从您的矩阵绘制,报告的行和列的ij

i j 
1 1 1 
2 2 2 
3 2 2