非负整数
设 n
为样本大小:
x <- rmultinom(n, V, rep.int(1 / G, G))
这是一个G x n
矩阵,其中每一列都是多项式样本,总和为V
。
通过将rep.int(1 / G, G)
传递给参数prob
,我假设每个组具有相等的“成功”概率。
正整数
正如Gregor所提到的,多项式样本可以包含0。如果不希望出现这种样本,则应该将其拒绝。因此,我们从截断的多项式分布中进行采样。
在如何在拒绝准则下从分布中生成目标数量的样本中,我建议采用“过度采样”方法来实现截断采样的“向量化”。简单来说,知道接受概率后,我们可以估计预期的试验次数M
来看到第一个“成功”(非零)。我们首先采样大约1.25 * M
个样本,然后这些样本中至少会有一个“成功”。我们随机返回其中一个作为输出。
以下函数实现了这个想法,以生成没有0的截断多项式样本。
positive_rmultinom <- function (n, V, prob) {
G <- length(prob)
if (G > V) stop("'G > V' causes 0 in a sample for sure!")
if (any(prob < 0)) stop("'prob' can not contain negative values!")
sum_prob <- sum(prob)
if (sum_prob != 1) prob <- prob / sum_prob
min_prob <- min(prob)
M <- round(1.25 * 1 / min_prob)
N <- n * M
x <- rmultinom(N, V, prob)
keep <- which(colSums(x == 0) == 0)
x[, sample(keep, n)]
}
现在让我们尝试一下。
V <- 76
prob <- c(53, 13, 9, 1)
直接使用rmultinom
绘制样本有时会导致其中一些样本为0:
sum(colSums(rmultinom(1000, V, prob) == 0) > 0)
但是使用positive_rmultinom
不存在这样的问题:
sum(colSums(positive_rmultinom(1000, V, prob) == 0) > 0)
combn
的确是个好点子。至于递增顺序,随机抽取列数不就解决了吗? - Brian Davism
有数值rep(1:(V - 1), G - 1)
(只有G - 1
列)。然后,在你展开网格之后,你可以添加最后一列以完成总和tmp = cbind(tmp, V - rowSums(tmp))
。你仍然需要取出添加的列小于1的行。但是,从expand.grid
结果中有效地删除一个维度将有助于使它能够扩展到比原来更大的规模。 - Gregor Thomas