生成总和为固定值的非负（或正）随机整数

Question

我想将正整数随机分配给G组，使它们的总和为V。

\n\n

例如，如果G = 3和V = 21，有效结果可能是(7, 7, 7)、(10, 6, 5)等。

\n\n

有没有一种简单的方法可以做到这一点？

\n\n

---

\n\n

编者注（来自\xe6\x9d\x8e\xe5\x93\xb2\xe6\xba\x90）：

\n\n

如果值不限于整数，则问题很简单，并且已在选择具有固定总和的 n 个数字中得到解决。

\n\n

对于整数，之前有一个问答：在R中生成N个随机整数，其总和为M，但它看起来更复杂并且很难理解。那边基于循环的解决方案也不令人满意。

\n

Answer 1

非负整数

设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) {
  ## 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