Generate N random integers that sum to M in R

I would like to generate N random positive integers that sum to M. I would like the random positive integers to be selected around a fairly normal distribution whose mean is M/N, with a small standard deviation (is it possible to set this as a constraint?).

Finally, how would you generalize the answer to generate N random positive numbers (not just integers)?

I found other relevant questions, but couldn't determine how to apply their answers to this context: 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

Solution 1:

Normalize.

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
}

For a continuous version, simply use:

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

Examples

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

Solution 2:

Just came up with an algorithm to generate N random numbers greater or equal to k whose sum is S, in an uniformly distributed manner. I hope it will be of use here!

First, generate N-1 random numbers between k and S - k(N-1), inclusive. Sort them in descending order. Then, for all x_i, with i <= N-2, apply x'_i = x_i - x_i+1 + k, and x'_N-1 = x_N-1 (use two buffers). The Nth number is just S minus the sum of all the obtained quantities. This has the advantage of giving the same probability for all the possible combinations. If you want positive integers, k = 0 (or maybe 1?). If you want reals, use the same method with a continuous RNG. If your numbers are to be integer, you may care about whether they can or can't be equal to k. Best wishes!

Explanation: by taking out one of the numbers, all the combinations of values which allow a valid Nth number form a simplex when represented in (N-1)-space, which lies at one vertex of a (N-1)-cube (the (N-1)-cube described by the random values range). After generating them, we have to map all points in the N-cube to points in the simplex. For that purpose, I have used one method of triangulation which involves all possible permutations of coordinates in descending order. By sorting the values, we are mapping all (N-1)! simplices to only one of them. We also have to translate and scale the numbers vector so that all coordinates lie in [0, 1], by subtracting k and dividing the result by S - kN. Let us name the new coordinates y_i.

Then we apply the transformation by multiplying the inverse matrix of the original basis, something like this:

    / 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 /

Which gives y'_i = y_i - y_i+1. When we rescale the coordinates, we get: 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, hence the above formula. This is applied to all elements except the last one.

Finally, we should take into account the distortion that this transformation introduces into the probability distribution. Actually, and please correct me if I'm wrong, the transformation applied to the first simplex to obtain the second should not alter the probability distribution. Here is the proof.

The probability increase at any point is the increase in the volume of a local region around that point as the size of the region tends to zero, divided by the total volume increase of the simplex. In this case, the two volumes are the same (just take the determinants of the basis vectors). The probability distribution will be the same if the linear increase of the region volume is always equal to 1. We can calculate it as the determinant of the transpose matrix of the derivative of a transformed vector V' = B^-1 V with respect to V, which, of course, is B^-1.

Calculation of this determinant is quite straightforward, and it gives 1, which means that the points are not distorted in any way that would make some of them more likely to appear than others.