生成 N 个随机整数，在 R 中总和为 M

Question

我想生成N 和为M 的随机正整数。我希望围绕一个相当正态的分布选择随机正整数，其平均值为M/N，具有较小的标准偏差（是否可以将其设置为约束？）。

最后，你将如何推广生成 N 个随机正数（不仅仅是整数）的答案？

我发现了其他相关问题，但无法确定如何将他们的答案应用于此上下文： https://stats.stackexchange.com/questions/59096/generate-three-random-numbers-that-sum-to-1-in-r

Generate 3 random number that sum to 1 in R

R - random approximate normal distribution of integers with predefined total

Answer 1

我想出了一个我认为更简单的解决方案。您首先从最小到最大范围生成随机整数，对它们进行计数，然后制作一个计数向量（包括零）。

请注意，即使最小值大于零，此解决方案也可能包含零。

希望这可以帮助未来的人解决这个问题:)

rand.vect.with.total <- function(min, max, total) {
  # generate random numbers
  x <- sample(min:max, total, replace=TRUE)
  # count numbers
  sum.x <- table(x)
  # convert count to index position
  out = vector()
  for (i in 1:length(min:max)) {
    out[i] <- sum.x[as.character(i)]
  }
  out[is.na(out)] <- 0
  return(out)
}

rand.vect.with.total(0, 3, 5)
# [1] 3 1 1 0

rand.vect.with.total(1, 5, 10)
#[1] 4 1 3 0 2

Answer 2

刚刚想出了一个算法，以均匀分布的方式生成 N 个大于或等于 k ​​且总和为 S 的随机数。希望对这里有用！

首先，在 k 和 S - k(N-1) 之间生成 N-1 个随机数，包括在内。按降序对它们进行排序。然后，对于所有 x_i，i i = x_i - x_{i+ 1} + k，并且 x'_N-1 = x_N-1（使用两个缓冲区）。第 N 个数字只是 S 减去所有获得量的总和。这具有为所有可能的组合提供相同概率的优点。如果你想要正整数，k = 0（或者可能是 1？）。如果您想要实数，请对连续 RNG 使用相同的方法。如果您的数字是整数，您可能会关心它们是否可以等于 k。祝你好运！

解释：通过取出其中一个数，所有允许有效第 N 个数的值的组合在 (N-1) 空间中表示时形成一个单纯形，该空间位于 (N-1) 的一个顶点处-cube（由随机值范围描述的 (N-1)-cube）。生成它们之后，我们必须将 N 立方体中的所有点映射到单纯形中的点。为此，我使用了一种三角测量方法，该方法涉及所有可能的坐标降序排列。通过对值进行排序，我们映射了所有 (N-1)！只简化为其中之一。我们还必须通过减去 k 并将结果除以 S - kN 来平移和缩放数字向量，以便所有坐标都位于 [0, 1] 中。让我们将新坐标命名为 y_i。

然后我们通过乘以原始基的逆矩阵来应用变换，如下所示：

    / 1  1  1 \            / 1 -1  0 \
B = | 0  1  1 |,    B^-1 = | 0  1 -1 |,    Y' = B^-1 Y
    \ 0  0  1 /            \ 0  0  1 /

这给出了 y'_i = y_i - y_i+1。当我们重新调整坐标时，我们得到： x'_i = y'_i(S - kN) + k = y_i(S - kN) - y_{i+ 1}(S - kN) + k = (x_i - k) - (x_i+1 - k) + k = x_i - x_i+1 + k，因此有上式。这适用于除最后一个以外的所有元素。

最后，我们应该考虑这种转换引入概率分布的失真。实际上，如果我错了，请纠正我，应用于第一个单纯形以获得第二个单纯形的变换不应该改变概率分布。这是证据。

任何点的概率增加是当区域大小趋于零时该点周围局部区域的体积增加除以单纯形的总体积增加。在这种情况下，两个体积是相同的（只取基向量的行列式）。如果区域体积的线性增加总是等于1，则概率分布将相同。我们可以将其计算为变换向量V' = B^{-1的导数的转置矩阵的行列式sup> V 相对于 V，当然是 B-1。}

此行列式的计算非常简单，它给出 1，这意味着这些点不会以任何方式扭曲，从而使其中一些点比其他点更容易出现。

Answer 3

标准化。

rand_vect <- function(N, M, sd = 1, pos.only = TRUE) {
  vec <- rnorm(N, M/N, sd)
  if (abs(sum(vec)) < 0.01) vec <- vec + 1
  vec <- round(vec / sum(vec) * M)
  deviation <- M - sum(vec)
  for (. in seq_len(abs(deviation))) {
    vec[i] <- vec[i <- sample(N, 1)] + sign(deviation)
  }
  if (pos.only) while (any(vec < 0)) {
    negs <- vec < 0
    pos  <- vec > 0
    vec[negs][i] <- vec[negs][i <- sample(sum(negs), 1)] + 1
    vec[pos][i]  <- vec[pos ][i <- sample(sum(pos ), 1)] - 1
  }
  vec
}

对于连续版本，只需使用：

rand_vect_cont <- function(N, M, sd = 1) {
  vec <- rnorm(N, M/N, sd)
  vec / sum(vec) * M
}

示例

rand_vect(3, 50)
# [1] 17 16 17

rand_vect(10, 10, pos.only = FALSE)
# [1]  0  2  3  2  0  0 -1  2  1  1

rand_vect(10, 5, pos.only = TRUE)
# [1] 0 0 0 0 2 0 0 1 2 0

rand_vect_cont(3, 10)
# [1] 2.832636 3.722558 3.444806

rand_vect(10, -1, pos.only = FALSE)
# [1] -1 -1  1 -2  2  1  1  0 -1 -1