Simulate a fair die with a biased die

Question

Given a biased $N$-sided die, how can a random number in the range $[1,N]$ be generated uniformly? The probability distribution of the die faces is not known, all that is known is that each face has a nonzero probability and that the probability distribution is the same on all throws (in particular, the throws are independent). This is the obvious generalization of Fair results with unfair die.

Putting this in computer science terms, we have an oracle representing the die rolls: $D : \mathbb{N} \to [1,N]$ such that $p_i = P(D(k)=i)$ is nonzero and independent of $k$. We're looking for a deterministic algorithm $A$ which is parametrized by $D$ (i.e. $A$ may make calls to $D$) such that $P(A()=i) = 1/N$. The algorithm must terminate with probability 1, i.e. the probability that $A$ makes more than $n$ calls to $D$ must converge to $0$ as $n\to\infty$.

For $N=2$ (simulate a fair coin from coin flips with a biased coin), there is a well-known algorithm:

• Repeat “flip twice” until the two throws come up with distinct outcomes ((heads, tails) or (tails, heads)). In other words, loop for $k = 0..\infty$ until $D(2k+1) \ne D(2k)$
• Return 0 if the last pair of flips was (heads, tails) and 1 if it was (tails, heads). In other words, return $D(2k)$ where $k$ is the index at which the loop was terminated.

A simplistic way to make an unbiased die from a biased one is to use the coin flip unbiasing method to build a fair coin, and build a fair die with rejection sampling, as in Unbiasing of sequences. But is this optimal (for generic values of the probability distribution)?

Specifically, my question is: what is an algorithm that requires the smallest expected number of calls to the oracle? If the set of reachable expected values is open, what is the lower bound and what is a class of algorithms that converges towards this lower bound?

In case different families of algorithms are optimal for different probability distributions, let's focus on almost-fair dice: I'm looking for an algorithm or a family of algorithms that's optimal for distributions such that $\forall i, \bigl|p_i - 1/N\bigr| \lt \epsilon$ for some $\epsilon \gt 0$.

Answer 1

The following paper answers a close variant of this question: The Efficient Construction of an Unbiased Random Sequence, Elias 1972.

The question there seems to be this: Given access to this biased independent source, output a sequence of random numbers in $[1,N]$ (note the difference from your question in which only one output symbol is requested). As the length of the desired output goes to infinity, the "efficiency" of the scheme in the paper (which seems like a natural generalization of von Neumann) goes to $1$, meaning, I believe, that an input with entropy $h$ is converted to an output of entropy approaching $h$.

The question seems much better behaved when phrased this way, rather than requesting a single output digit, because, for instance, if we draw $N$ samples and end up with an output with lots of information (for instance, all $N$ input symbols are distinct), then we can use all of that information to produce many output symbols, whereas with the question as phrased here, any information beyond that used to produce one output symbol goes to waste.

I believe that the scheme repeatedly takes $N$ draws, looks at the sequence, and maps it some outputs or the empty string. Perhaps there is a way to improve the scheme for your question by looking at prefixes and stopping if we have "enough" information to output a symbol? I don't know.

Answer 2

The method you describe for $N=2$ generalises. We use that all permutations of $[1..N]$ are equally likely even with a biased die (since the rolls are independent). Hence, we can keep rolling until we see such a permutation as the last $N$ rolls and output the last roll.

A general analysis is tricky; it is clear, however, that the expected number of rolls grows quickly in $N$ since the probability of seeing a permutation at any given step is small (and not independent of the steps before and after, hence tricky). It is greater than $0$ for fixed $N$, however, so the procedure terminates almost surely (i.e. with probability $1$).

For fixed $N$ we can construct a Markov chain over the set of Parikh-vectors that sum to $\leq N$, summarising the results of the last $N$ rolls, and determine the expected number of steps until we reach $(1,\dots,1)$ for the first time. This is sufficient since all permutations that share a Parikh-vector are equally likely; the chains and calculations are simpler this way.

Assume we are in state $v=(v_1,\dots,v_N)$ with $\sum_{i=1}^n v_i \leq N$. Then, the probability of gaining an element $i$ (i.e. the next roll is $i$) is always given by

$\qquad\displaystyle \operatorname{Pr}[\text{gain } i] = p_i$.

On the other hand, the propability of dropping an element $i$ from the history is given by

$\qquad\displaystyle \operatorname{Pr}_v[\text{drop } i] = \frac{v_i}{N}$

whenever $\sum_{i=1}^n v_i = N$ (and $0$ otherwise) precisely because all permutations with Parikh-vector $v$ are equally likely. These probabilities are independent (since the rolls are independent), so we can compute the transition probabilities as follows:

$\qquad\begin{align*} &\operatorname{Pr}[v \to (v_1,\dots,v_j + 1,\dots,v_N)] \\ &\qquad=\begin{cases} \operatorname{Pr}[\text{gain } j] &, \sum v < N \\ 0 &, \text{ else} \end{cases}\;,\\[3ex] &\operatorname{Pr}[v \to (v_1,\dots,v_i - 1, \dots v_j + 1,\dots,v_N)] \\ &\qquad=\begin{cases} 0 &, \sum v < N \lor v_i=0 \lor v_j=N \\ \operatorname{Pr}_v[\text{drop } i] \cdot \operatorname{Pr}[\text{gain } j] &, \text{ else} \end{cases}\;\text{ and}\\[3ex] &\operatorname{Pr}[v \to v] \\ &\qquad=\begin{cases} 0 &, \sum v < N \\ \sum_{v_i \neq 0} \operatorname{Pr}_v[\text{drop } i] \cdot \operatorname{Pr}[\text{gain } i] &, \text{ else} \end{cases}\;; \end{align*}$

all other transition probabilities are zero. The single absorbing state is $(1,\dots,1)$, the Parikh-vector of all permutations of $[1..N]$.

For $N=2$ the resulting Markov chain¹ is

^[source]

with expected number of steps until absorption

$\qquad\displaystyle \mathbb{E}\,\text{steps} = 2 p_0 p_1 \cdot 2 + \sum_{i\geq 3} (p_0^{i-1}p_1 + p_1^{i-1}p_0)\cdot i = \frac{1 - p_0 + p_0^2}{p_0-p_0^2}\;,$

using for simplification that $p_1=1-p_0$. If now, as suggested, $p_0 = \frac{1}{2} \pm \epsilon$ for some $\epsilon \in [0,\frac{1}{2})$, then

$\qquad\displaystyle \mathbb{E}\,\text{steps} = \frac{3 + 4\epsilon^2}{1 - 4\epsilon^2}$.

For $N \leq 6$ and uniform distributions (the best case) I have performed the calculations with computer algebra²; since the state space explodes quickly, larger values are hard to evaluate. The results (rounded upwards) are

^{Plots show $\mathbb{E}\,\text{steps}$ as a function of $N$; to the left a regular and to the right a logarithmic plot.}

The growth seems to be exponential but the values are too small to give good estimates.

As for stability against perturbations of the $p_i$ we can look at the situation for $N=3$:

^{Plot shows $\mathbb{E}\,\text{steps}$ as a function of $p_0$ and $p_1$; naturally, $p_2 = 1 - p_0 - p_1$.}

Assuming similar pictures for larger $N$ (kernel crashes computing symbolic results even for $N=4$), the expected number of steps seems to be quite stable for all but the most extreme choices (almost all or none mass at some $p_i$).

For comparison, simulating an $\epsilon$-biased coin (e.g. by assigning die results to $0$ and $1$ as evenly as possible), using this to simulate a fair coin and finally performing bit-wise rejection sampling requires at most

$\qquad\displaystyle 2\lceil \log N \rceil \cdot \frac{3 + 4\epsilon^2}{1 - 4\epsilon^2}$

die rolls in expectation -- you should probably stick with that.

---

1. Since the chain is absorbing in $(11)$ the edges hinted at in gray are never traversed and do not influence the calculations. I include them merely for completeness and illustrative purposes.
2. Implementation in Mathematica 10 (Notebook, Bare Source); sorry, it's what I know for these kinds of problems.

Answer 3

Just a quick comment regarding the case $N = 2$. Take some large number $m$, and sample $m$ throws of the die. If you got $k$ heads then you can extract $\log \binom{m}{k}$ bits. Assuming the die is $p$ biased, the average amount of information is $$ \sum_{k=0}^m p^k (1-p)^{m-k} \binom{m}{k} \log \binom{m}{k} \approx m h(p).$$ To get this estimate, use the fact that the binomial variable is concentrated around $k = pm$ together with the estimate $\log \binom{m}{k} \approx mh(k/m)$. As $m$ gets larger, we obtain the optimal rate of $h(p)$ per coin throw (this is optimal for information-theoretic reasons, for example the asymptotic equipartition property).

You can use the same method for general $N$, and you will probably get the same $H(\vec{p})$. These algorithms are only optimal in the limit, and there might be algorithms reaching the limit faster than these. In fact, I neglected to compute the speed of convergence - it might be an interesting exercise.

Answer 4

I would hazard the following answer.

The specific case of 2 you mentioned above is the specific case of expanding $(p + q)^2$ (where $p$ is prob of head and $q$ prob of tail) which gives you a term $2pq$ This means you can get $pq $ for one case and $qp$ for the other case. You will need to repeat sampling until you see either $pq$ or $qp$ (head-tail or tail-head) Using them as simulation, you will give equal probability.

When $N=3$ you have the expansion $(p + q + r)^3$ which gives you the term $pqr$. In this case, you do the same thing, sampling until you see all 3 outcomes $q$, $p$, $r$ in some order in 3 consecutive trials.

The same thing apply for the general case. Thinking carefully, I have to say the case of 2 is the best case where one can work things out in the expansion. When $N=3$ there are 6 different sequences for $pqr$ and there are many other terms in the expansion. I would feel quite uncomfortable with other terms where there are many more outcomes.

.

Extra:

This makes me think about the idea of simply sampling a lot to estimate the probability of each outcome of the dice. In this simplest case of one layer model with no hidden layer (a known model), we can work out a bound to conclude that the estimation converges quickly. In fact Chernoff bound shows us that the error goes down exponentially as sampling increases (linearly).

Now that a good estimation of the probabilities for each side of the dice is known, there are many options. One option is that we can do the expansion above again, but this time we can potentially use many other terms in the expansion that have the same value as $\prod_{i=1}^{i=n} p_i$ (or any term that you use as based sequence). This will be a bit more efficient because more terms in the expansion will be used. But I admit I don't know if this will result in the smallest number of calls to the oracle to have a guarantee on whatever preconditions (such as confidence parameter), if they are given.

Nevertheless, this approach is an answer to different flavor of the question. The question asks for guaranteed perfect unbiased-ness at the cost of potentially large sampling (though low prob). This approach only uses finite sampling with bound on confidence parameter. So I don't think this approach is appropriate to this question even though it is very interesting.