How to compute the distribution of sums when rolling 'N' dice with 'M' faces?

Question

I stumbled upon the following problem:

Given 'n' dice with 'm' faces with values 1 to m and a number 'x' what is the probability that the sum of the numbers on the 'm' dice is greater than or equal to 'x'? That is $m \le x \le n.m$ find $P(sum \ge x)$

Now, this 'almost approximates' a normal distribution. I'm sure it won't be a straight triangle but would have some sort of 'bell'-ish shape. Now the point is, "how" can I compute this without necessarily enumerating everything? Is there a generating function like thing I can use. Probably a recursive definition? Not sure if that's the way to go but I'm just at my wit's end. I was thinking of visualizing this distribution but that would still require to know the individual probabilities $P(sum = 6),...,P(sum = m.n)$

I've mostly encountered problems with 2-3 dice and it is somewhat easy to do or compute the sums. But how can one generalize this like in this problem?

Answer 1

For many purposes, you can use a normal approximation with continuity correction.

You can also use recursion for the probability that the total equals $y$, and add these from $y=x+1$ through $y=mn$.

There is also an exact formula for the probability that the sum equals $y$ involving a single summation, which gives you a double sum for the probability the sum is greater than $x$. Be careful that in the linked answer, the "dice" had values from $0$ through $m$ instead of from $1$ through $m$. You can convert these by subtracting $1$ from each die and $n$ from the total, and using $m-1$ in place of $m$. Also, as whuber pointed out, there can be numerical instability if you are not careful in the order in which you add the terms of an alternating sum.

Answer 2

Why don't you compute $P(sum > x)$ using a recursive formula with respect to the number of dices? Something like $P(sum_n > x) = P(sum_{n - 1} + outcome_n > x) = \sum_{1}^{m}p_i1_{1 \le i \le m}P(sum_{n-1}>x-i)$.

Answer 3

Here's an implementation of @ThePawn's answer in Julia. This computes the probability that they sum up to exactly n.

using Memoize # using Pkg; Pkg.add("Memoize")
@memoize function p_dice(dice, sides, n)
    # Returns the probability dice dice with side sides
    # sum up to n,
    # where side ∈ 1:sides
    if dice == 1
        if 1 <= n <= sides
            return 1/sides
        else
            return 0
        end
    end
    return sum(1/sides * p_dice(dice-1, sides, n-outcome) for outcome in 1:sides)
end

and one in python:

import functools
@functools.lru_cache
def p_dice(dice, sides, n):
    # Returns the probability dice dice with side sides
    # sum up to n,
    # where side in range(1, side+1)
    if dice == 1:
        if 1 <= n <= sides:
            return 1/sides # equal probability for each outcome
        else:
            return 0
    return sum(1/sides * p_dice(dice-1, sides, n-outcome) for outcome in range(1, sides+1)

Note that the memoization may very well be suboptimal for larger number of dice.

Answer 4

If the faces are uniform then the distribution will resemble a discrete form of the Irwin Hall distribution.

You could figure out the equations for the distribution, which will be some piecewise polynomial.

---

Possibly easier is to compute it by convolution (as given in another answer using a recursive formula) or by using the normal approximation.