我想将正整数随机分配给G
组,这样它们的总和为V
。
例如,如果 G = 3
和 V = 21
,则有效结果可能是 (7, 7, 7)
,(10, 6, 5)
等
有没有直接的方法来做到这一点?
编者注(来自 李哲源 ):
如果值不限于整数,问题很简单,已在 Choosing n numbers with fixed sum 中解决。 .
对于整数,之前有个问答:Generate N random integers that sum to M in R但它看起来更复杂并且很难理解。那边基于循环的解决方案也不令人满意。
最佳答案
非负整数
让 n
为样本大小:
x <- rmultinom(n, V, rep.int(1 / G, G))
是一个 G x n
矩阵,其中每一列都是一个 multinomial总和为 V
的示例。
通过将 rep.int(1/G, G)
传递给参数 prob
,我假设每个组都有相同的“成功”概率。
正整数
作为Gregor提到,多项式样本可以包含 0。如果此类样本不受欢迎,则应将其拒绝。因此,我们从截断的多项式分布中采样。
在How to generate target number of samples from a distribution under a rejection criterion我建议采用“过度采样”方法来实现截断采样的“矢量化”。简而言之,了解接受概率后,我们可以估算出预期的试验次数 M
才能看到第一个“成功”(非零)。我们首先采样说 1.25 * M
个样本,那么这些样本中至少会有一个“成功”。我们随机返回一个作为输出。
下面的函数实现了这个想法来生成不带 0 的截断多项式样本。
positive_rmultinom <- function (n, V, prob) {
## input validation
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!")
## normalization
sum_prob <- sum(prob)
if (sum_prob != 1) prob <- prob / sum_prob
## minimal probability
min_prob <- min(prob)
## expected number of trials to get a "success" on the group with min_prob
M <- round(1.25 * 1 / min_prob)
## sampling
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:
## number of samples that contain 0 in 1000 trials
sum(colSums(rmultinom(1000, V, prob) == 0) > 0)
#[1] 355 ## or some other value greater than 0
但是使用positive_rmultinom
就没有这个问题了:
## number of samples that contain 0 in 1000 trials
sum(colSums(positive_rmultinom(1000, V, prob) == 0) > 0)
#[1] 0
关于r - 生成总和为固定值的非负(或正)随机整数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/52559455/