# Algorithm for rolling an infinite-sided weighted die

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 any1 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 at2.

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.

• One thing I thought of is some variant of "generate some fraction 0 < x <= 1; choose 1/x", or 2/x or 3/x, but I'd like to see other answers with exponential decay. Jul 8, 2015 at 22:41
• Are you familiar with the exponential distribution? See en.wikipedia.org/wiki/Exponential_distribution Jul 8, 2015 at 22:42
• Please note that $60 + 25 + 12 + 4 + 1 + 0.2 + 0.04 > 100$, so the probabilities you have chosen are not possible in a discrete probability distribution. en.wikipedia.org/wiki/… Jul 9, 2015 at 4:59

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.

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$.

• This is how I usually like to think of it, but it doesn't quite work out as you present it - more precisely, one needs to find the solution in $n$ for $g(n)\leq k < g(n+1)$. It cannot be the case that $g$ has an inverse, given that $\mathbb N$ and $[0,1]$ have different cardinality. (But, nevertheless, the intuition that we're calculating an "inverse-like" (precisely, a left-inverse) of $g$ is useful) Jul 8, 2015 at 22:56
• Oops, i had that in my head when i was writing this but then my brain turned off and my hands took on the job Jul 8, 2015 at 22:57

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.

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$.

• this is a nice intuition, but computationally inadvisable Jul 9, 2015 at 12:18