Algorithm for rolling an infinite-sided weighted die

Question

If I wanted to have a die that rolled, for example:

| Roll | Prob (in %) |
|------|-------------|
| 1    | 60          |
| 2    | 25          |
| 3    | 12          |
| 4    | 4           |
| 5    | 1           |
| 6    | 0.2         |
| 7    | 0.04        |
| ...  | ...         |

(These numbers are off the top of my head, and don't necessarily follow a simple algorithm)

What kinds of algorithms can I programmatically describe such that I could roll any¹ number with exponentially decreasing odds?

I was thinking that if I just generated a table of probabilities using an exponential curve I wouldn't be able to guarantee my function resolves; If I started at the 1 roll and had a few failures in a row, I might never have a succeeding roll; and I can't start at the other end at work my way towards a roll of 1 with a probability of 100%, because it's infinitely long and there is no end to start at².

So how can I pick a number from this distribution? If possible, the distribution should exponentially or quadratically decay, but any distribution that could conceivably model something like:

1. The number of rooms in a house
2. The number of items in a treasure chest
3. The number of eyeballs on a mutant

A final but not necessary point of consideration would be the ease of which the algorithm can be modified to stretch/shorten the curve, to make low numbers less/more likely.

---

^{1: Barring computer precision; a roll of MAXINT or whatever would be the practical maximum
2: Again, in practice a computer could do it because there is a smallest float for the probability to be, but that would be horribly inefficient. The average case would take billions of rolls, whereas a good algorithm could probably do it in O(1) by just picking some fraction and resolving it.}

Answer 1

These probabilities are described by a Geometric Distribution. For that one, the probability of outcome $n$ is given by: $$P(X = n) = (1-p)^{n-1}p$$ where $p$ is the "exponential rate".

To sample from this distribution (or other ones for that matter) you can do the following. You choose from a uniform $[0,1]$ random variable and then look at the inverse value of that outcome through the CDF. The CDF for this distribution is $$P(X \leq n) = 1 - (1-p)^n$$ So here you would find the first $n$ s.t. that function is greater than the uniform outcome you chose.

Answer 2

If you have an $f: \mathbb{N} \rightarrow [0,1]$ function such that $n$ is a side and $f(n)$ the probability of getting that side, and if we can produce a function $g: \mathbb{N} \rightarrow [0,1]$ such that $g$ is invertible and:

$$ g(n) = \sum_{i=1}^n f(n) $$

Then you could just produce a random float $ k\in [0,1]$ and return $g^{-1}(k)$

Edit: As was pointed out, we can't expect $g^{-1}(k)$ to exist. However, maybe if we pretend for a minute that $g$ is defined at $\mathbb{R} \ge 0 $, and if we're lucky and $g'(x) \ge 0, \forall x \in \mathbb{R} \ge 0 $, then we can just do $\left\lfloor g^{-1}(k)\right\rfloor$.

Answer 3

The standard approach here is to use the CDF of your distribution, which is defined on $\mathbb{R}$ (regardless of the support of your distribution).

For any finite distribution, you generally use cumsum to calculate the CDF; for an infinite distribution, you'll need an explicit formula for $F:\mathbb{R}\rightarrow[0,1]$, where $F(x)=\mathbb{P}[X\leq x]$.

For a discrete distribution, the CDF will be a right-continuous step function, e.g., for the geometric distribution (image Wikipedia):

Once you have your CDF, draw a uniform normal variate $u$~$U[0,1]$. We need to invert $F$, which is done by picking the least $x$ s.t. $f(x)$ (the density) is defined and $F(x)\geq u$, i.e. we pick the $x$ that gets as close to $u$ without going under.

To see why, consider a simple coin flip, represented by a Bernoulli(1/4) variable. We want to go from $u\in[0,1]$ to something that shows up 0 a quarter of the time and 1 the other three quarters. Note that our cumulative distribution is $F(x)=\frac{1}{4}\mathbb{1}[x\geq 0] + \frac{3}{4}\mathbb{1}[x \geq 1]$, where $\mathbb{1}[\cdot]$ is an indicator function taking the value 1 if its argument is true, 0 otherwise.

By the rule I described, we'll pick $0 \Leftrightarrow u<\frac{1}{4}$ and $1 \Leftrightarrow u \geq \frac{1}{4}$; note that this means we'll pick 0 a quarter of the time, and 1 three quarters, which is what we wanted.

Put one more way, the probability of hitting any given value is given by the distance on the $y$ axis between the jumps in the CDF at that value. Distances on the $y$ axis correspond to probabilities for $U[0,1]$ variates, so to match the probability of a value, we should assign that value to $u$ whenever $u$ falls in that gap.

Answer 4

The inverse transform method mentioned by two of the other answers is one of the most elegant for generating random numbers over distributions. The exponential distribution has a very simple inverse CDF that makes numerical generation easy; you can apply the ceiling function to get whole-number outputs.

An alternative algorithm that you can use with real dice on a tabletop is simply "roll an n-sided die; if the number rolled is less than or equal to m then roll again, otherwise stop; the result is the number of rolls before stopping.", which obeys a geometric distribution. It gives a result of 1 with probability $1 - {m\over n}$, 2 with probability ${m \over n} \left(1 - {m \over n}\right)$, 3 with probability ${m \over n}^2 \left(1 - {m \over n}\right)$, etc. It has an expected (average) result of $n \over {n - m}$, which is the same as the expected number of rolls to get the result. The procedure terminates with probability 1, because for any probability $p$, no matter how small, there is a finite number of rolls $r$ such that the probability of not terminating after $r$ rolls is less than $p$.