Accurately simulating the lots of dice rolls without loops?

Question

OK so if your game rolls lots of dice you can just call a random number generator in a loop. But for any set of dice being rolled often enough you will get a distribution curve/histogram. So my question is there a nice simple calculation I can run that will give me a number that fits that distribution?

E.g. 2D6 - Score - % Probability

2 - 2.77%

3 - 5.55%

4 - 8.33%

5 - 11.11%

6 - 13.88%

7 - 16.66%

8 - 13.88%

9 - 11.11%

10 - 8.33%

11 - 5.55%

12 - 2.77%

So knowing the above you could roll a single d100 and work out an accurate 2D6 value. But once we start with 10D6, 50D6, 100D6, 1000D6 this could save a lot of processing time. So there must be a tutorial / method / algorithm that can do this fast? It is probably handy for stock markets, casinos, strategy games, dwarf fortress etc. What if you could simulate the outcomes of complete strategic battle that would take hours to play with a few calls to this function and some basic maths?

Answer 1

As I mentioned in my comment above, I recommend you profile this before overcomplicating your code. A quick for loop summing dice is a lot easier to understand and modify than complicated math formulae and table-building/searching. Always profile first to make sure you're solving the important problems. ;)

That said, there are two main ways to sample sophisticated probability distributions in one fell swoop:

---

1. Cumulative Probability Distributions

There's a neat trick to sample from continuous probability distributions by using only a single uniform random input. It has to do with the cumulative distribution, the function that answers "What is the probability of getting a value no greater than x?"

This function is non-decreasing, starting at 0 and rising to 1 over its domain. An example for the sum of two six-sided dice is shown below:

If your cumulative distribution function has a convenient-to-calculate inverse (or you can approximate it with piecewise functions like Bézier curves), you can use this to sample from the original probability function.

The inverse function handles parcelling up the domain between 0 and 1 into intervals mapped to each output of the original random process, with the catchment area of each matching its original probability. (This is true infinitescimally for continuous distributions. For discrete distributions like dice rolls we need to apply careful rounding)

Here's an example of using this to emulate 2d6:

int SimRoll2d6()
{
    // Get a random input in the half-open interval [0, 1).
    float t = Random.Range(0f, 1f);
    float v;

    // Piecewise inverse calculated by hand. ;)
    if(t <= 0.5f)
    {
         v = (1f + sqrt(1f + 288f * t)) * 0.5f;
    }
    else
    {
         v = (25f - sqrt(289f - 288f * t)) * 0.5f;
    }

    return floor(v + 1);
}

Compare this to:

int NaiveRollNd6(int n)
{
    int sum = 0;
    for(int i = 0; i < n; i++)
       sum += Random.Range(1, 7); // I'm used to Range never returning its max
    return sum;
}

See what I mean about the difference in code clarity and flexibility? The naive way might be naive with its loops, but it's short and simple, immediately obvious about what it does, and easy to scale to different die sizes and numbers. Making changes to the cumulative distribution code requires some non-trivial math, and it would be easy to break and cause unexpected results without any obvious mistakes. (Which I hope I haven't made above)

So, before you do away with a clear loop, make absolutely certain that it's really a performance problem worth this kind of sacrifice.

---

2. The Alias Method

The cumulative distribution method works well when you can express the inverse of the cumulative distribution function as a simple math expression, but that's not always easy or even possible. A reliable alternative for discrete distributions is something called the Alias Method.

This lets you sample from any arbitrary discrete probability distribution by using just two independent, uniformly distributed random inputs.

It works by taking a distribution like the one below on the left (don't worry that the areas/weights don't sum to 1, for the Alias Method we care about relative weight) and converting it to a table like the one on the right where:

• There is one column for each outcome.
• Each column is split into at most two parts, each associated with one of the original outcomes.
• The relative area/weight of each outcome is preserved.

(Diagram based on images from this excellent article on sampling methods)

In code, we represent this with two tables (or a table of objects with two properties) representing the probability of choosing the alternate outcome from each column, and the identity (or "alias") of that alternate outcome. Then we can sample from the distribution like so:

int SampleFromTables(float[] probabiltyTable, int[] aliasTable)
{
    int column = Random.Range(0, probabilityTable.Length);
    float p = Random.Range(0f, 1f);
    if(p < probabilityTable[column])
    {
        return column;
    }
    else
    {
        return aliasTable[column];
    }
}

This involves a bit of set-up:

1. Compute the relative probabilities of every possible outcome (so if you're rolling 1000d6, we need to compute the number of ways to get every sum from 1000 to 6000)

2. Build a pair of tables with an entry for each outcome. The full method goes beyond the scope of this answer, so I highly recommend referring to this explanation of the Alias Method algorithm.

3. Store those tables and refer back to them each time you need a new random die roll from this distribution.

This is a space-time tradeoff. The precomputation step is somewhat exhaustive, and we need to set aside memory proportionate to the number of outcomes we have (though even for 1000d6, we're talking single-digit kilobytes, so nothing to lose sleep over), but in exchange our sampling is constant-time no matter how complex our distribution might be.

---

I hope one or the other of those methods may be of some use (or that I've convinced you that the naive method's simplicity is worth the time it takes to loop) ;)

Answer 2

The answer unfortunately is that this method would not result in an increase in performance.

I believe that there may be some misunderstanding in the question of how a random number is generated. Take the example below [Java]:

Random r = new Random();
int n = 20;
int min = 1; //arbitrary
int max = 6; //arbitrary
for(int i = 0; i < n; i++){
    int randomNumber = (r.nextInt(max - min + 1) + min)); //silly maths
    System.out.println("Here's a random number: " + randomNumber);
}

This code will loop 20 times printing out random numbers between 1 and 6 (inclusive). When we talk about the performance of this code, there is some time taken to create the Random object (which involves creating an array of pseudo-random integers based on the computer's internal clock at the time it was created), and then 20 constant time look-ups on each nextInt() call. Since each 'roll' is a constant time operation, this makes rolling very cheap time-wise. Also notice that the range of min to max doesn't matter (in other words, it's just as easy for a computer to roll a d6 as it is for it to roll a d10000). Speaking in terms of time complexity, the performance of the solution is simply O(n) where n is the number of dice.

Alternatively, we could approximate any number of d6 rolls with a single d100 (or d10000 for that matter) roll. Using this method, we first need to compute s [number of faces to the dice] * n [number of dice] percentages before we roll (technically it's s * n - n + 1 percentages, and we should be able to divide that roughly in half since it's symmetrical; notice that in your example for simulating a 2d6 roll, you calculated 11 percentages and 6 were unique). After rolling we can use a binary search to figure out which range our roll fell into. In terms of time complexity, this solution evaluates to an O(s*n) solution, where s is the number of sides and n is the number of dice. As we can see, this is slower than the O(n) solution proposed in the previous paragraph.

Extrapolating from there, say you created both of these programs to simulate a 1000d20 roll. The first would simply roll 1,000 times. The second program would first need to determine 19,001 percentages (for the potential range of 1,000 to 20,000) before doing anything else. So unless you're on a weird system where memory look-ups are severally more expensive than floating-point operations, using a nextInt() call for each roll seems like the way to go.

Answer 3

If you want to store the dice combinations, the good news is that there is a solution, the bad is that our computers are somehow limited in regards to this kind of problems.

The good news:

There is a determistic approach of this problem:

1/ Compute all the combinations of your group of dice

2/ Determine the probability for each combination

3/ Search in this list for a result instead of throwing the dices

The bad news:

The number of combination with repetition is given by the following formulae

\$\Gamma_n^k={n+k-1 \choose k}=\frac{(n+k-1)!}{k!~(n-1)!}\$

(from french wikipedia):

That means that, for example, with 150 dices, you have 698'526'906 combinations. Let's assume you store the probability as a 32 bit float, you will need 2,6 GB of memory, and you still have to add memory requirement for the indexes...

In computing terms, the combination number can be computed by convolutions, which is handy but doesn't solve the memory constraints.

In conclusion, for a high number of dices, I would advise throwing the dices and observe the result rather than precomputing the probabilities associated with each combination.

Edit

However, as you are only interested in the sum of dices, you can store the probabilities with much less resources.

You can calculate precise probabilities for each sum of dice using convolution.

The general formula is \$F_i(m) = \sum_n F_1(n) F_{i-1}(m - n)\$

Then starting from 1/6 form each result with 1 dice, you can construct all the correct probabilities for any number of dice.

Here is a rough java code I wrote for illustration (not really optimized) :

public class DiceProba {

private float[][] probas;
private int currentCalc;

public int getCurrentCalc() {
    return currentCalc;
}

public float[][] getProbas() {
    return probas;
}

public void calcProb(int faces, int diceNr) {

    if (diceNr < 0) {
        currentCalc = 0;
        return;
    }

    // Initialize
    float baseProba = 1.0f / ((float) faces);
    probas = new float[diceNr][];
    probas[0] = new float[faces + 1];
    probas[0][0] = 0.0f;
    for (int i = 1; i <= faces; ++i)
        probas[0][i] = baseProba;

    for (int i = 1; i < diceNr; ++i) {

        int maxValue = (i + 1) * faces + 1;
        probas[i] = new float[maxValue];

        for (int j = 0; j < maxValue; ++j) {

            probas[i][j] = 0;
            for (int k = 0; k <= j; ++k) {
                probas[i][j] += probability(faces, k, 0) * probability(faces, j - k, i - 1);
            }

        }

    }

    currentCalc = diceNr;

}

private float probability(int faces, int number, int diceNr) {

    if (number < 0 || number > ((diceNr + 1) * faces))
        return 0.0f;

    return probas[diceNr][number];

}

}

Call calcProb() with the parameters you want and then access the proba table for results (first index: 0 for 1 dice, 1 for two dices...).

I checked it with 1'000D6 on my laptop, it took 10 seconds to compute all probabilities from 1 to 1'000 dices and all possible sums of dices.

With precomputing and efficient storage, you could have quick answers for a high number of dices.

Hope it helps.