在R中生成N个随机整数,使它们的和为M

22
我希望生成N个随机正整数,使它们的和为M。我希望这些随机正整数是在一个均值为M/N、标准偏差较小的正态分布内选择的(是否可以将其设置为约束条件?)。
最后,如何将答案概括以生成N个随机正数(不仅仅是整数)?
我找到了其他相关的问题,但无法确定如何将它们的答案应用到这个上下文中: https://stats.stackexchange.com/questions/59096/generate-three-random-numbers-that-sum-to-1-in-r 在R中生成三个总和为1的随机数 R-随机近似正态分布的整数,预定义总数

我还没有阅读那些文章,但它们听起来确实相关。 - Strawberry
我认为我没有理解这个问题并完全欣赏下面的解决方案。这里有一个更整洁的问答:生成总和为固定值的非负(或正)随机整数。希望对本帖的读者有所帮助。 - Zheyuan Li
3个回答

25

标准化。

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

能否从均匀分布中生成N个随机整数,使它们的和为M? - Nell
不同的问题。看起来你可能会在Rhelp上找到一个几乎相同的问题的答案。(在这里和Rhelp上都不建议交叉发布。) - IRTFM
1
rand_vect_cont不能保证输出仅为正数。例如,请参见set.seed(1984) rand_vect_cont(3, 10) - Lubo

2
刚刚想出一个算法,可以以均匀分布的方式生成N个大于等于k且总和为S的随机数。希望这对您有用!
首先,生成N-1个介于k和S-k(N-1)之间(包括边界)的随机数。将它们按降序排序。然后,对于所有xi,其中i≤N-2,应用x'i = xi - xi+1 + k,并且x'N-1 = xN-1(使用两个缓冲区)。第N个数字就是S减去所有获得量的总和。这样做的优点是为所有可能的组合提供相同的概率。如果您想要正整数,则k = 0(或者也许是1?)。如果您想要实数,则使用具有连续RNG的相同方法。如果您的数字必须是整数,则可能会关心它们是否能够等于k。祝一切顺利!
解释:通过取出其中一个数字,允许有效的第N个数字的所有值的组合在表示为(N-1)空间中时形成一个单形体,该单形体位于(N-1)个立方体中的一个顶点处(随机值范围描述的(N-1)个立方体)。生成它们后,我们必须将N-立方体中的所有点映射到单形体中的点。为此,我使用了一种三角剖分方法,其中包括坐标在降序中的所有可能的排列。通过对值进行排序,我们将所有(N-1)!个单形体映射到其中的一个。我们还必须将数值向量进行平移和缩放,使得所有坐标都在[0, 1]内,通过减去k并将结果除以S-kN 来实现。让我们将新坐标命名为yi
然后,我们通过乘以原始基础的逆矩阵来应用变换,类似于这样:
    / 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 = yi - yi+1。当我们重新调整坐标时,得到以下公式:

x'i = y'i(S - kN) + k = yi(S - kN) - yi+1(S - kN) + k = (xi - k) - (xi+1 - k) + k = xi - xi+1 + k,因此上述公式适用于除了最后一个元素之外的所有元素。

最后,我们应该考虑这种转换引入概率分布的扭曲。实际上,如果我说错了,请纠正我,用于获得第二个简单形式的第一个简单形式的转换不应该改变概率分布。以下是证明。

在任何点处的概率增加是当该区域的大小趋近于零时,该点周围局部区域的体积增加与简单形式总体积增加的比值。在这种情况下,两个体积相同(只需取基向量的行列式)。如果区域体积的线性增长始终等于1,则概率分布将保持不变。我们可以将其计算为一个转换向量V' = B-1 V的导数的转置矩阵的行列式,这当然是B-1

计算这个行列式非常简单,它等于1,这意味着点没有扭曲成使一些点比其他点更有可能出现的任何方式。


2
我建议你重新发明一些接近对称狄利克雷分布的东西。去查一下吧,维基百科上有一些漂亮的图形。 - IRTFM

0

我找到了一个我认为更简单的解决方案。首先从最小值到最大值范围内生成随机整数,计算它们的数量,然后制作一个包括零的计数向量。

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

希望这能帮助未来遇到这个问题的 R 语言用户 :)

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

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接