How to generate a non-integer amount of consecutive Bernoulli successes?

Question

Given:

1. A coin with unknown bias $p$ (Head).
2. A strictly positive real $a > 0$.

Problem:

Generate a random Bernoulli variate with bias $p^{a}$.

Does anyone know how to do this? For instance, when $a$ is a positive integer, then one can flip the coin $a$ times and see whether all the outcomes were Heads: if they are then issue '0', otherwise issue '1'. The difficulty lies in the fact that $a$ is not necessarily an integer. Also, if I knew the bias $p$, I could just build another coin with the desired bias.

Answer 1

We can solve this via a couple of "tricks" and a little math.

Here is the basic algorithm:

1. Generate a Geometric random variable with probability of success $p$.
2. The outcome of this random variable determines a fixed known value $f_n \in [0,1]$.
3. Generate a $\mathrm{Ber}(f_n)$ random variable using fair coin flips generated from blockwise paired flips of our $\mathrm{Ber}(p)$ coin.
4. The resulting outcome will be $\mathrm{Ber}(p^a)$ for any $a \in (0,1)$, which is all we need.

To make things more digestible, we'll break things into pieces.

Piece 1: Without loss of generality assume that $0 < a < 1$.

If $a \geq 1$, then, we can write $p^a = p^n p^b$ for some positive integer $n$ and some $0 \leq b < 1$. But, for any two independent Bernoulli's, we have $$\renewcommand{\Pr}{\mathbb P} \Pr(X_1 = X_2 = 1) = p_1 p_2 \>. $$ We can generate a $p^n$ Bernoulli from our coin in the obvious way. Hence, we need only concern ourselves with generating the $\mathrm{Ber}(p^a)$ when $a \in (0,1)$.

Piece 2: Know how to generate an arbitrary $\mathrm{Ber}(q)$ from fair coin flips.

There is a standard way to do this. Expand $q = 0.q_1 q_2 q_3 \ldots$ in its binary expansion and then use our fair coin flips to "match" the digits of $q$. The first match determines whether we declare a success ("heads") or failure ("tails"). If $q_n = 1$ and our coin flip is heads, declare heads, if $q_n = 0$ and our coin flip is tails, declare tails. Otherwise, consider the subsequent digit against a new coin flip.

Piece 3: Know how to generate a fair coin flip from unfair ones with unknown bias.

This is done, assuming $p \in (0,1)$, by flipping the coin in pairs. If we get $HT$, declare a heads; if we get $TH$, declare a tails, and otherwise repeat the experiment until one of the two aforementioned outcomes occurs. They are equally probable, so must have probability $1/2$.

Piece 4: Some math. (Taylor to the rescue.)

By expanding $h(p) = p^a$ around $p_0 = 1$, Taylor's theorem asserts that $$ p^a = 1 - a(1-p) - \frac{a(1-a)}{2!} (1-p)^2 - \frac{a(1-a)(2-a)}{3!} (1-p)^3 \cdots \>. $$ Note that because $0 < a < 1$, each term after the first is negative, so we have $$ p^a = 1 - \sum_{n=1}^\infty b_n (1-p)^n \>, $$ where $0 \leq b_n \leq 1$ are known a priori. Hence $$ 1 - p^a = \sum_{n=1}^{\infty} b_n (1-p)^n = \sum_{n=1}^\infty b_n \Pr(G \geq n) = \sum_{n=1}^\infty f_n \Pr(G = n) = \mathbb E f(G), $$ where $G \sim \mathrm{Geom}(p)$, $f_0 = 0$ and $f_n = \sum_{k=1}^n b_k$ for $n \geq 1$.

And, we already know how to use our coin to generate a Geometric random variable with probability of success $p$.

Piece 5: A Monte Carlo trick.

Let $X$ be a discrete random variable taking values in $[0,1]$ with $\Pr(X = x_n) = p_n$. Let $U \mid X \sim \mathrm{Ber}(X)$. Then $$ \Pr(U = 1) = \sum_n x_n p_n. $$

But, taking $p_n = p(1-p)^n$ and $x_n = f_n$, we see now how to generate a $\mathrm{Ber}(1-p^a)$ random variable and this is equivalent to generating a $\mathrm{Ber}(p^a)$ one.

Answer 2

Is the following answer silly?

If $X_1,\dots,X_n$ are independent $\mathrm{Ber}(p)$ and $Y_n$ has distribution $\mathrm{Ber}\left(\left(\sum_{i=1}^n X_i/n \right)^a\right)$, then $Y_n$ will be approximately distributed as $\mathrm{Ber}(p^a)$, when $n\to\infty$.

Hence, if you don't know $p$, but you can toss this coin a lot of times, it is possible to sample (approximately) from a $\mathrm{Ber}(p^a)$ random variable.

Example Rcode:

n <- 1000000
p <- 1/3 # works for any 0 <= p <= 1
a <- 4
x <- rbinom(n, 1, p)
y <- rbinom(n, 1, mean(x)^a)
cat("p^a =", p^a, "\n")
cat("est =", mean(y))

Results:

p^a = 0.01234568 
est = 0.012291 

Answer 3

I posted the following exposition of this question and cardinal's answer to the General Discussion forum of the current Analytic Combinatorics class on Coursera, "Application of power series to constructing a random variable." I'm posting a copy here as community wiki to make this publicly and more permanently available.

---

There was an interesting question and answer on stat.stackexchange.com related to power series: "How to generate a non-integer amount of consecutive Bernoulli successes?" I'll paraphrase the question and the answer by cardinal.

Suppose we have a possibly unfair coin which is heads with probability $p$, and a positive real number $\alpha$. How can we construct an event whose probability is $p^\alpha$?

If $\alpha$ were a positive integer, we could just flip the coin $\alpha$ times and let the event be that all tosses were heads. However, if $\alpha$ is not an integer, say $1/2$, then this doesn't make sense, but we can use this idea to reduce to the case that $0 \lt \alpha \lt 1$. If we want to construct an event whose probability is $p^{3.5}$, we take the intersection of independent events whose probabilities are $p^3$ and $p^{0.5}$.

One thing we can do is construct an event with any known probability $p' \in [0,1]$. To do this, we can construct a stream of fair bits by repeatedly flipping the coin twice, reading $HT$ as $1$ and $TH$ as $0$, and ignoring $HH$ and $TT$. We compare this stream with the binary expansion of $p' = 0.a_1a_2a_3..._2$. The event that the first disagreement is where $a_i=1$ has probability $p'$. We don't know $p^\alpha$, so we can't use this directly, but it will be a useful tool.

The main idea is that we would like to use the power series for $p^\alpha = (1-q)^\alpha = 1 - \alpha q - \frac{\alpha(1-\alpha)}{2} q^2 - \frac{\alpha (1-\alpha)(2-\alpha)}{3!}q^3 -...$ where $p=1-q$. We can construct events whose probabilities are $q^n$ by flipping the coin $n$ times and seeing if they are all tails, and we can produce an event with probability $p' q^n$ by comparing the binary digits of $p'$ with a fair bit stream as above and checking whether $n$ tosses are all tails.

Construct a geometric random variable $G$ with parameter $p$. This is the number of tails before the first head in an infinite sequence of coin tosses. $P(G=n) = (1-p)^np = q^n p$. (Some people use a definition which differs by $1$.)

Given a sequence $t_0, t_1, t_2, ...$, we can produce $t_G$: Flip the coin until the first head, and if there are $G$ tails before the first head, take the element of the sequence of index $G$. If each $t_n \in [0,1]$, we can compare $t_G$ with a uniform random variable in $[0,1]$ (constructed as above) to get an event with probability $E[t_G] = \sum_n t_n P(G=n) = \sum_n t_n q^n p $.

This is almost what we need. We would like to eliminate that $p$ to use the power series for $p^\alpha$ in $q$.

$$1 = p + qp + q^2p + q^3p + ...$$

$$q^n = q^np + q^{n+1}p + q^{n+2}p + ...$$

$$\begin{eqnarray} \sum_n s_n q^n & = & \sum_n s_n (q^n p + q^{n+1}p + q^{n+2}p + ...) \newline & = & \sum_n (s_0 + s_1 + ... + s_n) q^n p \end{eqnarray}$$

Consider $1-p^\alpha = \alpha q + \frac{\alpha(1-\alpha)}{2} q^2 + ... $. Let $t_n$ be the sum of the coefficients of $q$ through $q^n$. Then $1-p^\alpha = \sum_n t_n q^n p$. Each $t_n\in [0,1]$ since the coefficients are positive and sum to $1-0^\alpha = 1$, so we can construct an event with probability $1-p^\alpha$ by comparing a fair bit stream with the binary expansion of $t_G$. The complement has probability $p^\alpha$ as required.

---

Again, the argument is due to cardinal.

Answer 4

The very complete answer by cardinal and subsequent contributions inspired the following remark/variant.

Let PZ stand "Probability of Zero" and $q:=1-p$. If $X_n$ is an iid Bernoulli sequence with PZ $q$, then $M_n := \max(X_1,\,X_2,\,\dots, X_n)$ is a Bernoulli r.v. with PZ $q^n$. Now making $n$ random i.e., replacing it by an integer rv $N \geq 1$ leads to Bernoulli rv $M_N$ with $$ \mathrm{Pr}\{M_N =0\} = \sum_{n=1}^\infty \mathrm{Pr}\{M_N =0 \,\vert\, N =n\} \mathrm{Pr}\{N =n\} = \sum_{n=1}^\infty \mathrm{Pr}\{N =n\} \, q^n. $$ So if $0 < a < 1$ and if we take $\mathrm{Pr}\{N =n\} =b_n$ from cardinal's answer, we find $\mathrm{Pr}\{M_N =0\} = 1- p^a$ and $1-M_N$ is $\mathrm{Ber}(p^a)$ as wanted. This is indeed possible since the coefficients $b_n$ satisfy $b_n \geqslant 0$ and they sum to $1$.

The discrete distribution of $N$ depends only on $a$ with $0 < a < 1$, recall $$ \mathrm{Pr}\{N =n\} = \frac{a}{n}\,\prod_{k=1}^{n-1}\left(1 - a/k\right) \qquad (n \geq 1). $$ It has interesting features. It turns out to have an infinite expectation and an heavy tail behaviour $n \,b_n \sim c/n^a$ with $c = -1/\Gamma(-a) >0$.

Though $M_N$ is the maximum of $N$ rvs, its determination needs a number of $X_k$ which is $\leq N$ since the result is known as soon as one $X_k$ is $1$. The number of computed $X_k$ is geometrically distributed.