Can someone tell me how to simulate $mathrm{Bernoulli}left({aover b}right)$, where $a,bin mathbb{N}$, using a coin toss (as many times as you require) with $P(H)=p$ ?

I was thinking of using rejection sampling, but could not nail it down.

#### Best Answer

**Because there are uncountably many solutions, let's find an efficient one.**

The idea behind this one starts with a standard way to implement a Bernoulli variable: compare a uniform random variable $U$ to the parameter $a/b$. When $U lt a/b$, return $1$; otherwise, return $0$.

**We can use the $p$-coin as a uniform random number generator**. To generate a number $U$ uniformly within any interval $[x, y)$, flip the coin. When it's heads, recursively generate a uniform value $X$ in the first $p$ part of the interval; when it's tails, recursively generate $X$ from the last $1-p$ part of the interval. At some point the target interval will become so small that it doesn't really matter how you pick a number from it: that's how the recursion gets started. It's obvious this procedure generates uniform variates (up to any desired precision), as is easily proven by induction.

**This idea is not efficient, but it leads to an efficient method.** Since at each stage you are going to draw a number from some given interval $[x,y)$, why not first check whether you need to draw it at all? If your target value lies outside this interval, *you already know the result of the comparison between the random value and the target.* Thus, this algorithm tends to terminate rapidly. (This could be construed as the *rejection sampling* procedure requested in the question.)

**We can optimize this algorithm further.** At any stage, we actually have *two* coins we can use: by relabeling our coin we can make it into one that is heads with chance $1-p$. Therefore, as a precomputation *we may recursively choose whichever relabeling leads to the lower expected number of flips needed for termination.* (This calculation can be an expensive step.)

For example, it's inefficient to use a coin with $p=0.9$ to emulate a Bernoulli$(0.01)$ variable directly: it takes almost ten flips on average. But if we use a $p=1-0.0=0.1$ coin, then in just two flips we will sure to be done and the expected number of flips is just $1.2$.

**Here are the details.**

Partition any given half-open interval $I = [x, y)$ into the intervals

$$[x,y) = [x, x + (y-x)p) cup [x + (y-x)p, y) = s(I,H) cup s(I,T).$$

This defines the two transformations $s(*,H)$ and $s(*,T)$ which operate on half-open intervals.

As a matter of terminology, if $I$ is any set of real numbers let the expression

$$t lt I$$

mean that $t$ is a lower bound for $I$: $t lt x$ for all $x in I$. Similarly, $t gt I$ means $t$ is an upper bound for $I$.

Write $a/b = t$. (In fact, it will make no difference if $t$ is real instead of rational; we only require that $0 le t le 1$.)

**Here is the algorithm to produce a variate $Z$ with the desired Bernoulli parameter:**

Set $n=0$ and $I_n = I_0 = [0,1)$.

While $(tin I_{n})$ {Toss the coin to produce $X_{n+1}$. Set $I_{n+1} = S(I_n, X_{n+1}).$ Increment $n$.}

If $t gt I_{n+1}$ then set $Z=1$. Otherwise, set $Z=0$.

### Implementation

To illustrate, here is an `R`

implementation of the alorithm as the function `draw`

. Its arguments are the target value $t$ and the interval $[x,y)$, initially $[0,1)$. It uses the auxiliary function `s`

implementing $s$. Although it does not need to, also it tracks the number of coin tosses. It returns the random variable, the count of tosses, and the last interval it inspected.

`s <- function(x, ab, p) { d <- diff(ab) * p if (x == 1) c(ab[1], ab[1] + d) else c(ab[1] + d, ab[2]) } draw <- function(target, p) { between <- function(z, ab) prod(z - ab) <= 0 ab <- c(0,1) n <- 0 while(between(target, ab)) { n <- n+1; ab <- s(runif(1) < p, ab, p) } return(c(target > ab[2], n, ab)) } `

As an example of its use and test of its accuracy, take the case $t=1/100$ and $p=0.9$. Let's draw $10,000$ values using the algorithm, report on the mean (and its standard error), and indicate the average number of flips used.

`target <- 0.01 p <- 0.9 set.seed(17) sim <- replicate(1e4, draw(target, p)) (m <- mean(sim[1, ])) # The mean (m - target) / (sd(sim[1, ]) / sqrt(ncol(sim))) # A Z-score to compare to `target` mean(sim[2, ]) # Average number of flips `

In this simulation $0.0095$ of the flips were heads. Although lower than the target of $0.01$, the Z-score of $-0.5154$ is not significant: this deviation can be attributed to chance. The average number of flips was $9.886$–a little less than ten. If we had used the $1-p$ coin, the mean would have been $0.0094$–still not significantly different than the target, but only $1.177$ flips would have been needed on average.

### Similar Posts:

- Solved – Given Bernoulli probability, how to draw a Bernoulli from a uniform distribution
- Solved – Given Bernoulli probability, how to draw a Bernoulli from a uniform distribution
- Solved – Given Bernoulli probability, how to draw a Bernoulli from a uniform distribution
- Solved – Given Bernoulli probability, how to draw a Bernoulli from a uniform distribution
- Solved – Generating discrete uniform from coin flips