57
$\begingroup$

This is a question I found on Glassdoor: How does one generate 7 integers with equal probability using a coin that has a $\mathbb{Pr}(\text{Head}) = p\in(0,1)$?

Basically, you have a coin that may or may not be fair, and this is the only random-number generating process you have, so come up with random number generator that outputs integers from 1 to 7 where the probability of getting each of these integers is 1/7.

Efficiency of the data-generates process matters.

$\endgroup$
  • $\begingroup$ Related, non-dupe: stats.stackexchange.com/questions/256563/… $\endgroup$ – Sycorax Jul 5 '18 at 20:29
  • 12
    $\begingroup$ There are myriad ways to accomplish this. A more interesting version of the question asks for the best method in some well-defined sense. A natural sense of best would be least expected number of flips per integer generated. Another interesting version is to describe all possible solutions (that rely on independent flips of the coin and nothing more). $\endgroup$ – whuber Jul 5 '18 at 21:16
  • 1
    $\begingroup$ @whuber good suggestion, I've edited the question to reflect your comment. $\endgroup$ – Amazonian Jul 6 '18 at 0:52
  • $\begingroup$ <<<Basically, you have a coin that may or may not be fair, and this is the only random-number generating process you have>>> Does this mean, using the coin in any different method than flipping it and checking for head vs. tail is "forbidden", since it would be another random-number generating process? $\endgroup$ – TinglTanglBob Jul 6 '18 at 10:46
  • 9
    $\begingroup$ Mod 7 of the year on the coin. $\endgroup$ – Nat Jul 6 '18 at 19:25

11 Answers 11

55
$\begingroup$

Flip the coin twice. If it lands HH or TT, ignore it and flip it twice again.

Now, the coin has equal probability of coming up HT or TH. If it comes up HT, call this H1. If it comes up TH, call this T1.

Keep obtaining H1 or T1 until you have three in a row. These three results give you a number based on the table below:

H1 H1 H1 -> 1
H1 H1 T1 -> 2
H1 T1 H1 -> 3
H1 T1 T1 -> 4
T1 H1 H1 -> 5
T1 H1 T1 -> 6
T1 T1 H1 -> 7
T1 T1 T1 -> [Throw out all results so far and repeat]

I argue that this would work perfectly fine, although you would have a lot of wasted throws in the process!

$\endgroup$
  • 4
    $\begingroup$ The only constraint is that the probability of heads is "p". Note that p could be $0$, or $1$. So this isn't guaranteed to work, but Sycorax's (or Sephan's) would work, even in those cases. $\endgroup$ – gung Jul 5 '18 at 20:34
  • 8
    $\begingroup$ @gung: I'm not sure I'd work for a coin with two heads, or two tails. $\endgroup$ – Stephan Kolassa Jul 5 '18 at 21:01
  • 6
    $\begingroup$ With probability $1$ it will finish in finite time. $\endgroup$ – clark Jul 5 '18 at 23:15
  • 18
    $\begingroup$ This is called Von-Neuman whitening. $\endgroup$ – DonFusili Jul 6 '18 at 8:37
  • 4
    $\begingroup$ You can iterate the von Neumann extractor to more fully extract entropy from the sequence. Gather all the HH and TT pairs, consider that a sequence, apply the von Neumann extractor to that, etc.. $\endgroup$ – A Simmons Jul 6 '18 at 21:00
46
$\begingroup$

Assume that $p \in (0,1)$.

Step 1:. Toss the coin 5 times.

If the outcome is

$(H, H, H, T, T)$, return $1$ and stop.

$(H, H, T, T, H)$, return $2$ and stop.

$(H, T, T, H, H)$, return $3$ and stop.

$(T, T, H, H, H)$, return $4$ and stop.

$(T, H, H, H, T)$, return $5$ and stop.

$(H, H, T, H, T)$, return $6$ and stop.

$(H, T, H, T, H)$, return $7$ and stop.

Step 2:. If the outcome is none of the above, repeat Step 1.

Note that regardless of the value of $p \in (0,1)$, each of the seven outcomes listed above has probability $q = p^3(1-p)^2$, and the expected number of coin tosses is $\displaystyle \frac{5}{7q}$. The tosser doesn't need to know the value of $p$ (except that $p\neq 0$ and $p\neq 1$); it is guaranteed that the seven integers are equally likely to be returned by the experiment when it terminates (and it is guaranteed to end with probability $1$).

$\endgroup$
  • 6
    $\begingroup$ Can we reduce the expected number of tosses for this by allowing either the sequence specified here OR that sequence with each flip inverted. eg: For 1, either (H,H,H,T,T) or (T,T,T,H,H)? $\endgroup$ – moreON Jul 6 '18 at 6:40
  • 5
    $\begingroup$ You could add the complement as well. If the outcome is (H,H,H,T,T) or (T,T,T,H,H), return 1 and stop, etc. In that case the probability for each outcome is $q=p^3(1-p)^2+p^2(1-p)^3$. $\endgroup$ – Martijn Weterings Jul 6 '18 at 8:45
  • 2
    $\begingroup$ Wouldn't it be possible to just add another coin-flip if the result is not withing any (H,H,H,T,T) arrangement? With the additional coin-toss you'd need another mapping of (H,H,H,T,T,T) and (H,H,T,T,T,T) and each xT(7-x)H combination that can be arranged in 7 or more different orders to the numbers 1 to 7. Instead of retossing all 5 coins this would add only 1 additional toss, but im not sure if it works :D $\endgroup$ – TinglTanglBob Jul 6 '18 at 13:14
  • 5
    $\begingroup$ Maybe it might be the best thing to instantly flip the coin 7 times, since than it is guaranteed, that you will get a random number out of it (only exception therefore is, that the coin lands heads up or tail up all 7 tries). So with 7 tosses you could end up with 1 up to 6 heads (i exclude the 0 and 7 option here since its pointless). If one head there are 7 different arrangements of (H,T,T,T,T,T,T) possible; If 2 heads its 21; if 3 heads its 35; if 4 35 heads; if 5 21 heads; if 6 7 heads; Each can be perfectly mapped to number 1-7 without any combination lost. $\endgroup$ – TinglTanglBob Jul 6 '18 at 13:47
  • 2
    $\begingroup$ @TinglTanglBob This is essentially Martijn Weterings answer ;-) $\endgroup$ – M.Herzkamp Jul 6 '18 at 14:01
20
$\begingroup$

Divide a box into seven equal-area regions, each labeled with an integer. Throw the coin into the box in such a way that it has equal probability of landing in each region.

This is only half in jest -- it's essentially the same procedure as estimating $\pi$ using a physical Monte Carlo procedure, like dropping rice grains onto paper with a circle drawn on it.

This is one of the only answer that works for the case of $p = 1$ or $p=0$.

$\endgroup$
  • 2
    $\begingroup$ Can you suggest a way to divide the box into seven equal-area regions, to minimize bias from flipping, bouncing off walls etc? Seven sectors of angle 360/7 ? $\endgroup$ – smci Jul 5 '18 at 23:35
  • 1
    $\begingroup$ @smci That's why I stipulate that you must throw the coin so that there is uniform probability of it landing in each square. If bouncing off of a wall influences that probability, then you have to account for it in your throw. $\endgroup$ – Sycorax Jul 6 '18 at 0:00
  • 17
    $\begingroup$ Yes I know that, and I'm pointing out to you that simply saying "throw it in an unbiased way" without defining exactly how to achieve that, is not really a full answer... in which case the flipping-H/T based methods are superior. $\endgroup$ – smci Jul 6 '18 at 0:03
  • 1
    $\begingroup$ "Equal areas, with throws that have equal probability of landing in each region" might be difficult to establish in practice. In practice it might be easier to mark out a large number of "practice runs" and then subdivide the landing area into equiprobable spaces empirically (eg with 700 throws, rule off a line that cuts off the furthest 100 throws, then another for the next 100 and so on). A simplification of this to generate a single random bit would be to throw the coin twice - if the first throw goes further, the bit is 0, and if the second throw goes further it's a 1 $\endgroup$ – Silverfish Jul 6 '18 at 9:54
  • 4
    $\begingroup$ There's a nice answer by @TheScienceBoy that sadly is deleted with an interesting alternative to this - effectively, using the coin as a spinner, and marking 7 sections along its circumference - that retains much of the spirit of this answer but may be more physically straightforward to carry out! $\endgroup$ – Silverfish Jul 6 '18 at 9:56
19
$\begingroup$

Generalizing the case described by Dilip Sarwate

Some of the methods described in the other answers use a scheme in which you throw a sequence of $n$ coins in a 'turn' and depending on the result you choose a number between 1 or 7 or discard the turn and throw again.

The trick is to find in the expansion of possibilities a multiple of 7 outcomes with the same probability $p^k(1-p)^{n-k}$ and match those against each other.

Because the total number of outcomes is not a multiple of 7, we have a few outcomes that we can not assign to a number, and have some probability that we need to discard the outcomes and start over.


The case of using 7 coin flips per turn

Intuitively we could say that rolling the dice seven times would be very interesting. Since we only need to throw out 2 out of the $2^7$ possibilities. Namely, the 7 times heads and 0 times heads.

For all other $2^7-2$ possibilities there is always a multiple of 7 cases with the same number of heads. Namely 7 cases with 1 heads, 21 cases with 2 heads, 35 cases with 3 heads, 35 cases with 4 heads, 21 cases with 5 heads, and 7 cases with 6 heads.

So if you compute the number (discarding 0 heads and 7 heads) $$X = \sum_{k=1}^{7} (k-1) \cdot C_k $$

with $C_k$ Bernoulli distributed variables (value 0 or 1), then X modulo 7 is a uniform variable with seven possible results.


Comparing different number of coin flips per turn

The question remains what the optimal number of rolls per turn would be. Rolling more dices per turn cost you more, but you reduce the probability to have to roll again.

The image below shows a manual computations for the first few numbers of coin flips per turn. (possibly there might be an analytical solution, but I believe it is safe to say that a system with 7 coin flips provides the best method regarding the expectation value for the necessary number of coin flips)

expected number of coin flips

# plot an empty canvas
plot(-100,-100,
     xlab="flips per turn",
     ylab="E(total flips)",
     ylim=c(7,400),xlim=c(0,20),log="y")
title("expectation value for total number of coin flips
(number of turns times flips per turn)")

# loop 1
# different values p from fair to very unfair 
# since this is symmetric only from 0 to 0.5 is necessary 

# loop 2
# different values for number of flips per turn
# we can only use a multiple of 7 to assign 
#   so the modulus will have to be discarded
#   from this we can calculate the probability that the turn succeeds
#   the expected number of flips is 
#       the flips per turn 
#             divided by 
#       the probability for the turn to succeed 

for (p in c(0.5,0.2,0.1,0.05)) {
  Ecoins <- rep(0,16)
  for (dr in (5:20)){
    Pdiscards = 0
    for (i in c(0:dr)) { 
      Pdiscards = Pdiscards + p^(i)*(1-p)^(dr-i) * (choose(dr,i) %% 7)
    }
    Ecoins[dr-4] = dr/(1-Pdiscards)
  }
  lines(5:20, Ecoins)
  points(5:20, Ecoins, pch=21, col="black", bg="white", cex=0.5)
  text(5, Ecoins[1], paste0("p = ",p), pos=2)
}

Using an early stopping rule

note: the calculations below, for the expectation value of number of flips, are for a fair coin $p=0.5$, it would become a mess to do this for different $p$, but the principle remains the same (although different book-keeping of the cases is needed)

We should be able to choose the cases (instead of the formula for $X$) such that we might be able to stop earlier.

  • With 5 coin flips we have for the six possible different unordered sets of heads and tails:

    1+5+10+10+5+1 ordered sets

    And we can use the groups with ten cases (that is the group with 2 heads or the group with 2 tails) to choose (with equal probability) a number. This occurs in 14 out of 2^5=32 cases. This leaves us with:

    1+5+3+3+5+1 ordered sets

  • With an extra (6-th) coin flip we have for the seven possible different unordered sets of heads and tails:

    1+6+8+6+8+6+1 ordered sets

    And we can use the groups with eight cases (that is the group with 3 heads or the group with 3 tails) to choose (with equal probability) a number. This occurs in 14 out of 2*(2^5-14)=36 cases. This leaves us with:

    1+6+1+6+1+6+1 ordered sets

  • With another (7-th) extra coin flip we have for the eight possible different unordered sets of heads and tails:

    1+7+7+7+7+7+7+1 ordered sets

    And we can use the groups with seven cases (all except the all tails and all heads cases) to choose (with equal probability) a number. This occurs in 42 out of 44 cases. This leaves us with:

    1+0+0+0+0+0+0+1 ordered sets

    (we could continue this but only in the 49-th step does this give us an advantage)

So the probability to select a number

  • at 5 flips is $\frac{14}{32} = \frac{7}{16}$
  • at 6 flips is $\frac{9}{16}\frac{14}{36} = \frac{7}{32}$
  • at 7 flips is $\frac{11}{32}\frac{42}{44} = \frac{231}{704}$
  • not in 7 flips is $1-\frac{7}{16}-\frac{7}{32}-\frac{231}{704} = \frac{2}{2^7}$

This makes the expectation value for the number of flips in one turn, conditional that there is success and p=0.5:

$$5 \cdot \frac{7}{16}+ 6 \cdot \frac{7}{32} + 7 \cdot \frac{231}{704} = 5.796875 $$

The expectation value for the total number of flips (until there is a success), conditional that p=0.5, becomes:

$$\left(5 \cdot \frac{7}{16}+ 6 \cdot \frac{7}{32} + 7 \cdot \frac{231}{704}\right) \frac{2^7}{2^7-2} = \frac{53}{9} = 5.88889 $$


The answer by NcAdams uses a variation of this stopping-rule strategy (each time come up with two new coin flips) but is not optimally selecting out all the flips.

The answer by Clid might be similar as well although there might be an uneven selection rule that each two coin flips a number might be chosen but not necessarily with equal probability (a discrepancy which is being repaired during later coin flips)


Comparison with other methods

Other methods using a similar principle are the one by NcAdams and AdamO.

The principle is: A decision for a number between 1 and 7 is made after a certain number of heads and tails. After an $x$ number of flips, for each decision that leads to a number $i$ there is a similar, equally probable, decision that leads to a number $j$ (the same number of heads and tails but just in a different order). Some series of heads and tails can lead to a decision to start over.

For such type of methods the one placed here is the most efficient because it makes the decisions as early as possible (as soon as there is a possibility for 7 equal probability sequences of heads and tails, after the $x$-th flip, we can use them to make a decision on a number and we do not need to flip further if we encounter one of those cases).

This is demonstrated by the image and simulation below:

comparison

#### mathematical part #####
set.seed(1)


#plotting this method
p <- seq(0.001,0.999,0.001)
tot <- (5*7*(p^2*(1-p)^3+p^3*(1-p)^2)+
       6*7*(p^2*(1-p)^4+p^4*(1-p)^2)+
       7*7*(p^1*(1-p)^6+p^2*(1-p)^5+p^3*(1-p)^4+p^4*(1-p)^3+p^5*(1-p)^2+p^6*(1-p)^1)+
        7*1*(0+p^7+(1-p)^7) )/
             (1-p^7-(1-p)^7)
plot(p,tot,type="l",log="y",
     xlab="p",
     ylab="expactation value number of flips"
     )

#plotting method by AdamO
tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
          (-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
          (21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
          35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
          1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
           *p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
           5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
           15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot,col=2,lty=2)

#plotting method by NcAdam
lines(p,3*8/7/(p*(1-p)),col=3,lty=2)

legend(0.2,500,
       c("this method calculation","AdamO","NcAdams","this method simulation"),
       lty=c(1,2,2,0),pch=c(NA,NA,NA,1),col=c(1,2,3,1))


##### simulation part ######

#creating decision table
mat<-matrix(as.numeric(intToBits(c(0:(2^5-1)))),2^5,byrow=1)[,c(1:12)]
colnames(mat) <- c("b1","b2","b3","b4","b5","b6","b7","sum5","sum6","sum7","decision","exit")

# first 5 rolls
mat[,8] <- sapply(c(1:2^5), FUN = function(x) {sum(mat[x,1:5])})

mat[which((mat[,8]==2)&(mat[,11]==0))[1:7],12] = rep(5,7) # we can stop for 7 cases with 2 heads
mat[which((mat[,8]==2)&(mat[,11]==0))[1:7],11] = c(1:7)   
mat[which((mat[,8]==3)&(mat[,11]==0))[1:7],12] = rep(5,7) # we can stop for 7 cases with 3 heads
mat[which((mat[,8]==3)&(mat[,11]==0))[1:7],11] = c(1:7)    

# extra 6th roll
mat <- rbind(mat,mat)
mat[c(33:64),6] <- rep(1,32)
mat[,9] <- sapply(c(1:2^6), FUN = function(x) {sum(mat[x,1:6])})

mat[which((mat[,9]==2)&(mat[,11]==0))[1:7],12] = rep(6,7) # we can stop for 7 cases with 2 heads
mat[which((mat[,9]==2)&(mat[,11]==0))[1:7],11] = c(1:7)   
mat[which((mat[,9]==4)&(mat[,11]==0))[1:7],12] = rep(6,7) # we can stop for 7 cases with 4 heads
mat[which((mat[,9]==4)&(mat[,11]==0))[1:7],11] = c(1:7)    

# extra 7th roll
mat <- rbind(mat,mat)
mat[c(65:128),7] <- rep(1,64)
mat[,10] <- sapply(c(1:2^7), FUN = function(x) {sum(mat[x,1:7])})

for (i in 1:6) {
  mat[which((mat[,10]==i)&(mat[,11]==0))[1:7],12] = rep(7,7) # we can stop for 7 cases with i heads
  mat[which((mat[,10]==i)&(mat[,11]==0))[1:7],11] = c(1:7)   
}


mat[1,12] = 7           # when we did not have succes we still need to count the 7 coin tosses
mat[2^7,12] = 7


draws = rep(0,100)
num = rep(0,100)
# plotting simulation
for (p in seq(0.05,0.95,0.05)) {
  n <- rep(0,1000)
  for (i in 1:1000) {
    coinflips <- rbinom(7,1,p)  # draw seven numbers
    I <- mat[,1:7]-matrix(rep(coinflips,2^7),2^7,byrow=1) == rep(0,7)                      # compare with the table
    Imatch = I[,1]*I[,2]*I[,3]*I[,4]*I[,5]*I[,6]*I[,7]        # compare with the table 
      draws[i] <- mat[which(Imatch==1),11]                 # result which number
      num[i]   <- mat[which(Imatch==1),12]                 # result how long it took
  }
  Nturn <- mean(num)                   #how many flips we made
  Sturn <- (1000-sum(draws==0))/1000   #how many numbers we got (relatively)
  points(p,Nturn/Sturn)
}

another image which is scaled by $p*(1-p)$ for better comparison:

comparison with scaled expectation values

zoom in comparing methods described in this post and comments

comparison methods described here

the 'conditional skipping of the 7-th step' is a slight improvement which can be made on the early stopping rule. In this case you select not groups with equal probabilities after the 6-th flips. You have 6 groups with equal probabilities, and 1 groups with a slightly different probability (for this last group you need to flip one more extra time when you have 6 heads or tails and because you discard the 7 heads or 7 tails, you will end up with the same probability after all)

$\endgroup$
  • $\begingroup$ I was just about starting to make the calculations for the n=7 case, because I had a feeling it might be better than n=1. Have my upvote, sir! $\endgroup$ – M.Herzkamp Jul 6 '18 at 11:39
  • $\begingroup$ @M.Herzkamp a small improvement is still possible. One number $C_k$ (one coin flip) is not necessary for the calculation of $X$, because it has coefficient 0. This number is only necessary to determine the case of all heads or all tails, and it can be omitted when we already know that we have a mixed case. $\endgroup$ – Martijn Weterings Jul 6 '18 at 13:08
  • $\begingroup$ So the improvement brings it down to a bit more than 6 coin flips as expectation value for the number of required coin flips. It would be vary though to proof that this is the optimal solution. The scheme created by Clid diverges a bit allowing to choose a number at a particular number of coin flips but not with equal probability (at least not for that particular step, it will be corrected later). $\endgroup$ – Martijn Weterings Jul 6 '18 at 16:38
  • $\begingroup$ But if you're deciding whether to flip the sixth coin based on the first five results, are the probabilities of each set, conditioned on you got to six flips, the same as for the other sets? $\endgroup$ – Acccumulation Jul 6 '18 at 17:56
  • $\begingroup$ @Acccumulation you could draw it as a binary tree with 7 levels. We will only be selecting among nodes if there are 7 with equal probability. It is like you chop of some branches earlier (at level 5 or 6). If you like then you could continue untill 7 steps, instead of earlier, but for these particular cases the 6th and 7th coin flip do not make a difference. $\endgroup$ – Martijn Weterings Jul 6 '18 at 18:07
7
$\begingroup$

EDIT: based on others' feedback.

Here's an interesting thought:

set the list of {1,2,3,4,5,6,7}. Throw the coin for each element in the list sequentially. If it lands head side up for a particular element, remove the number from the list. If all the numbers from a particular iteration of the list are removed, repeat the sampling. Do so until only one number remains.

drop.one <- function(x, p) {
  drop <- runif(length(x)) < p
  if (all(drop))
    return(x)
  return(x[!drop])
}

sample.recur <- function(x, p) {
  if (length(x) > 1)
    return(sample.recur(drop.one(x, p), p))
  return(x)
}

# x <- c(1:7,7:1)
x <- 1:7
p <- 0.01

out <- replicate(1e5, sample.recur(x, p))

round(prop.table(table(out)), 2)

gives me an approximately uniform distribution

> round(prop.table(table(out)), 2)
out
   1    2    3    4    5    6    7 
0.14 0.14 0.15 0.14 0.14 0.14 0.14 

It's interesting to note (if I haven't made a dire mistake) that this produces a different result than generating $N$ binomial outcomes as the sum of 13 tosses of the coin (counting 0 heads as an outcome) and mapping the {0,1,2,...,12,13} index onto the earlier list of {1,2,3,...,3,2,1}. I don't quite know how to prove that my method works.


Evaluation of expectation value for number of coin throws

The expectation value for the number coin throws can be calculated using the transition matrix below (answering the question when we start with $x$ non-eliminated numbers then what is the probability to get to $y$ non-eliminated numbers)

$$M = \begin{bmatrix} q^7 & 0 & 0 & 0 &0 & 0 & 1 & 1\\ 7p^1q^6 & q^6 & 0 & 0 & 0 & 0 & 0 & 0 \\ 21p^2q^5 & 6p^1q^5 & q^5 & 0 & 0 & 0 & 0 & 0\\ 35 p^3q^4 & 15 p^2q^4 & 5q^4 & q^4 & 0 & 0 & 0 & 0\\ 35 p^4q^3 & 20 p^3q^3 & 10 p^2q^3 & 4 p^1q^3 & q^3 & 0 & 0 & 0\\ 21p^5q^2 & 15 p^4q^2 & 10 p^3q^2 & 6 p^2q^2 & 3 p^1q^2 & q^2 & 0 & 0\\ 7p^6q^1 & 6 p^5q^1 & 5 p^4q^1 & 4 p^3q^1 & 3p^2q^1 & 2p^1q^1 & 0 & 0\\ p^7 & p^6 & p^5 & p^4 & p^3 & p^2 & 0 & 0 \end{bmatrix}$$

The eigenvector associated with the eigenvalue 1 (which can be found by solving $(M-I)v=0$) depicts how much time is relatively spend in what state. Then 7th state is how often you will be able to draw a number from 1 to 7. The other states tell how many coin flips it costs.

Below is the image which compares with the answer from NcAdams which has expectation value for coin throws being $E(n) = \frac{24}{7}p(1-p)$

comparison of expectation value for coin flips

Remarkable is that the method performs better roughly for $p>2/3$. But also the performance is non-symmetric. A symmetric and better overall performance could be made when a probabilistic switching rule would be made which changes the decision rule from tails to heads when heads happens to be improbable.

Solution found with wxMaxima

M: matrix(
 [(1-p)^7,        0,          0,0,0,0,1,1], 
 [7* p*(1-p)^6,   (1-p)^6,        0,0,0,0,0,0], 
 [21*p^2*(1-p)^5, 6*p*(1-p)^5,    (1-p)^5,0,0,0,0,0], 
 [35*p^3*(1-p)^4, 15*p^2*(1-p)^4, 5*p*(1-p)^4,(1-p)^4,0,0,0,0], 
 [35*p^4*(1-p)^3, 20*p^3*(1-p)^3, 10*p^2*(1-p)^3,4*p*(1-p)^3,(1-p)^3,0,0,0], 
 [21*p^5*(1-p)^2, 15*p^4*(1-p)^2, 10*p^3*(1-p)^2,6*p^2*(1-p)^2,3*p*(1-p)^2,(1-p)^2,0,0], 
 [7* p^6*(1-p)^1, 6*p^5*(1-p),    5*p^4*(1-p),4*p^3*(1-p),3*p^2*(1-p),2*(1-p)*p,0,0], 
 [p^7,        p^6,        p^5,p^4,p^3,p^2,0,0]
);
z: nullspace(M-diagmatrix(8,1));
x : apply (addcol, args (z));
t : [7,6,5,4,3,2,0,0];
plot2d(t.x/x[7],[p,0,1],logy);

Calculations in R

# plotting empty canvas
plot(-100,-100,
     xlab="p",
     ylab="E(total flips)",
     ylim=c(10,1000),xlim=c(0,1),log="y")

# plotting simulation
for (p in seq(0.1,0.9,0.05)) {

  n <- rep(0,10000)
  for (i in 1:10000) {
    success  = 0
    tests = c(1,1,1,1,1,1,1)     # start with seven numbers in the set
    count = 0
    while(success==0) {
      for (j in 1:7)  {
        if (tests[j]==1) {
          count = count + 1
          if  (rbinom(1,1,p) == 1) {
            tests[j] <- 0        # elliminate number when we draw heads
          }
        }
      }
      if (sum(tests)==1) {
        n[i] = count
        success = 1              # end     when 1 is left over
      }
      if (sum(tests)==0) {
        tests = c(1,1,1,1,1,1,1) # restart when 0 are left over
      }
    }
  }
  points(p,mean(n))
}

# plotting formula
p <- seq(0.001,0.999,0.001)

tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
    (-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
    (21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
    35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
    1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
   *p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
 5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
  15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot)

#plotting comparison with alternative method
lines(p,3*8/7/(p*(1-p)),lty=2)

legend(0.2,500,
       c("simulation","calculation","comparison"),
       lty=c(0,1,2),pch=c(1,NA,NA))
$\endgroup$
  • 1
    $\begingroup$ Clever idea (+1). Intuitively, it should work, since the symmetry would appear to nullify bias towards any particular number. Still, I'd love to see a proof. $\endgroup$ – Ben Jul 6 '18 at 2:36
  • 6
    $\begingroup$ This idea is really nice, but with a high probability for head (knock out that number) i think the last number in the row hast the best chances to "survive" since all the others infront get kicked out very likely before on run 1? Maybe it could be changed by not sequentially throwing the coin but parallel for all the numbers in x? Runtime of the script might increase i guess :) $\endgroup$ – TinglTanglBob Jul 6 '18 at 9:30
  • 2
    $\begingroup$ I agree with @TinglTanglBob - when I set p <- 0.99 I get the output 0.89 0.02 0.02 0.02 0.02 0.02 0.02 $\endgroup$ – Silverfish Jul 6 '18 at 9:48
  • 6
    $\begingroup$ Wouldn't performing the elimination in 'rounds' fix the bias problem? Start with 7 numbers. Toss the coin for every remaining number, and eliminate the ones that throw head. If all the remaining numbers are eliminated in a round, scratch the results of that round and try again. I don't know how to prove it, but intuitively the order of the numbers no-longer matters in whether they are the 'winner' $\endgroup$ – Phil Jul 6 '18 at 14:07
  • 1
    $\begingroup$ @TinglTanglBob thanks for pointing that out. I played with the script a bit more and my proposed solution is: use the side with less probability to be the elimination criterion. Setting $p = 0.01$ strangely fixes the problem. I echo the desire for one better than I to offer some proof for/against this "solution". $\endgroup$ – AdamO Jul 6 '18 at 14:41
5
$\begingroup$

The question is a bit ambiguous, is it asking "generate a random integer equal or less than 7 with equal probability", or is it asking "generate 7 random integers with equal probability?" - but what is the space of integers?!?

I'll assume it's the former, but the same logic I'm applying can be extended to the latter case too, once that problem is cleared up.

With a biased coin, you can produce a fair coin by following the following procedure: https://en.wikipedia.org/wiki/Fair_coin#Fair_results_from_a_biased_coin

A number 7 or less can be written in binary as three {0,1} digits. So all one needs to do is follow the above procedure three times, and convert the binary number produced back to decimal.

$\endgroup$
  • 1
    $\begingroup$ comparing my answer with @NcAdams, it's clear I'm including 0 as a possible desirable outcome! $\endgroup$ – Cam.Davidson.Pilon Jul 5 '18 at 20:25
  • $\begingroup$ I don't get how your answer is different. If you include {0,0,0} -> 1 then what does {1,1,1} map to? There are 8 possibilities. $\endgroup$ – AdamO Jul 5 '18 at 22:01
  • 1
    $\begingroup$ 000 maps to 0, hence my comment on including 0 as a possible value. This was posted before the OP edit and almost simultaneously as NcAdams. $\endgroup$ – Cam.Davidson.Pilon Jul 5 '18 at 22:47
  • $\begingroup$ Half of this answer is a comment, and the actual content of the answer is link-only. Please spell out your actual answer rather than just linking to it. $\endgroup$ – Cubic Jul 10 '18 at 10:33
3
$\begingroup$

A solution that never wastes flips, which helps a lot for very-biased coins.

The disadvantage of this algorithm (as written, at least) is that it's using arbitrary-precision arithmetic. Practically, you probably want to use this until integer overflow, and only then throw it away and start over.

Also, you need to know what the bias is ... which you might not, say, if it is temperature-dependent like most physical phenomena.


Assuming the chance of heads is, say, 30%.

  • Start with the range [1, 8).
  • Flip your coin. If heads, use the left 30%, so your new range is [1, 3.1). Else, use the right 70%, so your new range is [3.1, 8).
  • Repeat until the entire range has the same integer part.

Full code:

#!/usr/bin/env python3
from fractions import Fraction
from collections import Counter
from random import randrange


BIAS = Fraction(3, 10)
STAT_COUNT = 100000


calls = 0
def biased_rand():
    global calls
    calls += 1
    return randrange(BIAS.denominator) < BIAS.numerator


def can_generate_multiple(start, stop):
    if stop.denominator == 1:
        # half-open range
        stop = stop.numerator - 1
    else:
        stop = int(stop)
    start = int(start)
    return start != stop


def unbiased_rand(start, stop):
    if start < 0:
        # negative numbers round wrong
        return start + unbiased_rand(0, stop - start)
    assert isinstance(start, int) and start >= 0
    assert isinstance(stop, int) and stop >= start
    start = Fraction(start)
    stop = Fraction(stop)
    while can_generate_multiple(start, stop):
        if biased_rand():
            old_diff = stop - start
            diff = old_diff * BIAS
            stop = start + diff
        else:
            old_diff = stop - start
            diff = old_diff * (1 - BIAS)
            start = stop - diff
    return int(start)


def stats(f, *args, **kwargs):
    c = Counter()
    for _ in range(STAT_COUNT):
        c[f(*args, **kwargs)] += 1

    print('stats for %s:' % f.__qualname__)
    for k, v in sorted(c.items()):
        percent = v * 100 / STAT_COUNT
        print('  %s: %f%%' % (k, percent))


def main():
    #stats(biased_rand)
    stats(unbiased_rand, 1, 7+1)
    print('used %f calls at bias %s' % (calls/STAT_COUNT, BIAS))


if __name__ == '__main__':
    main()
$\endgroup$
  • 3
    $\begingroup$ You can improve this process. Start the from the interval [0,1]; you can see your sequence of coin tosses as an infinite process that generates a single random number in $[0,1]$, uniformly. This random number has a unique base-7 expansion (apart from a probability-0 set where you have two expansions, one ending in $0000\dots$ and one in $6666\dots$). When the interval is narrow enough that you can determine its $k$th (base-7) digit uniquely, you output it. It's a while since I read about that, but if I remember correctly, this generates digits at the optimal rate by Shannon's entropy bound. $\endgroup$ – Federico Poloni Jul 7 '18 at 13:19
  • $\begingroup$ That's the same for a single output, right? Just better for multiple? I'd write it as diff *= 7 I think ... in fact, there's no particular need to use the same base for each attempt. $\endgroup$ – o11c Jul 7 '18 at 15:14
  • $\begingroup$ Yes, it's the same if you want a single output; it only improves efficiency if you want multiple ones. $\endgroup$ – Federico Poloni Jul 7 '18 at 15:51
  • $\begingroup$ +1 This is in the spirit of the arithmetic coding procedure.And if $p$ is given as a rational number, this could be easily implemented exactly, with (unbounded) integer arithmetic. But, of course, ir requires to know $p$. $\endgroup$ – leonbloy Jul 7 '18 at 23:09
  • $\begingroup$ This absolutely does waste flips if you want a single output. For a fair coin, the standard technique (roll the coin three times, and repeat if you get TTT) gives an expected number of 24/7 = 3+3/7 rolls. If you use this arithmetic-coding-style technique, you roll at least four times unless you get HHH or TTT, which gives you an expected number of more than 15/4 = 3+3/4 rolls. $\endgroup$ – Peter Shor Jul 8 '18 at 17:46
3
$\begingroup$

As mentioned in earlier comments, this puzzle relates to John von Neumann's 1951 paper "Various Techniques Used in Connection With Random Digits" published in the research journal of the National Bureau of Standards:

enter image description here

There is a wider literature about such problems that goes under the name of Bernoulli factory problems, that is, given a coin with tail probability $p$, how to simulate a coin with tail probability $f(p)$. If feasible, since some functions $f$ cannot be used as for instance $f(p)=\min\{1,2p\}$. Nacu and Peres (2005) study fast algorithms for solving [solvable] Bernoulli factory problems where fast means exponential decay of the tail of the distribution of the number $N$ of trials.

$\endgroup$
2
$\begingroup$

This also only works for $p \neq 1$ and $p \neq 0$.

We first turn the (possibly) unfair coin into a fair coin using the process from NcAdams answer:

Flip the coin twice. If it lands HH or TT, ignore it and flip it twice again.

Now, the coin has equal probability of coming up HT or TH. If it comes up HT, call this H1. If it comes up TH, call this T1.

Now we use the fair coin to generate a real number between $0$ and $1$ in binary. Let H1$= 1$ and T1 $= 0$. Start with the string 0., flip the coin and append the resulting digit to the at the end of the string. Repeat with the new string. For example, the sequence H1 H1 T1 would give you the number $0.110$.

$1/7$ is a repeating decimal, and with the right-hand side being in base 2 we have that:

$1/7 = 0.001 001 001 ...$

$2/7 = 0.010 010 010 ...$

$3/7 = 0.011 011 011 ...$

$4/7 = 0.100 100 100 ...$

$5/7 = 0.101 101 101 ...$

$6/7 = 0.110 110 110 ...$

Keep flipping the fair coin to generate the decimal digits until the digits of your sequence does not match one of the above sequences, then chose the number $n$ such that your generated number is less than $n/7$ and greater than $(n-1)/7$. Since each generated number is equally likely we have chosen a number between $1$ and $7$ with equal probability.

$\endgroup$
  • 1
    $\begingroup$ Ignoring that this method requires a lot of flips, it is interesting how the separation of the bits into 7 equally likely parts are generated. These parts are really equally likely? For instance the $\leq 1/7$ number turns up when you flip 0.000 or when your not flipping of 0.001001001 ends up on the wrong side, which is 1/8+1/16 of the cases, or do I see this wrong? $\endgroup$ – Martijn Weterings Jul 6 '18 at 12:12
  • 1
    $\begingroup$ Correction on the probability to flip $\leq 1/7$ that last number is not 1/16. So you flip 000 or 001000 or 001001000, etc which has probability $\sum_{k=1}^{\infty} (1/8)^k = \frac{1}{7}$ $\endgroup$ – Martijn Weterings Jul 6 '18 at 12:20
  • 1
    $\begingroup$ So the probability for termination within $n$ bits is $T(n) = \frac{2^n-6}{2^n}$ (for $n \geq 3$ otherwise it is zero) and then the expected value for the necessary number of bits is $$3+\sum_{n=3}^{\infty} 1-T(n) = 4.5$$ this is not an improvement in comparison to NcAdams solution which needs $\frac{8}{7}*3 \approx 3.42$ bits. (the main culprit is that after 3 bits this method has only a decision for 25% of the cases) $\endgroup$ – Martijn Weterings Jul 6 '18 at 12:41
2
$\begingroup$

Inspired by AdamO's answer, here is a Python solution that avoids bias:

def roll(p, n):
    remaining = range(1,n+1)
    flips = 0
    while len(remaining) > 1:
        round_winners = [c for c in remaining if random.choices(['H','T'], [p, 1.0-p]) == ['H']]
        flips += len(remaining)
        if len(round_winners) > 0:
            remaining = round_winners
        p = 1.0 - p
    return remaining[0], flips

There are two main changes here: The main one is that if all the number are discarded in a round, repeat the round. Also I flip the choice of whether heads or tails means discard every time. This reduces the number of flips needed in cases where p is close to 0 or 1 by ~70% when p=0.999

$\endgroup$
  • 2
    $\begingroup$ "I flip the choice of whether heads or tails means discard every time. This reduces the number of flips needed in cases where p is close to 0 or 1 by ~70% when p=0.999" - smart thinking! $\endgroup$ – Silverfish Jul 6 '18 at 16:23
  • 1
    $\begingroup$ Alternating heads or tails is definitely an improvement over always discarding heads -- but perhaps it would be better if, after flipping a coin for each remaining option, if they are all the same we repeat reflipping all of them, otherwise if there are at least as many heads as tails we eliminate the remaining options corresponding to heads, otherwise we eliminate the remaining options corresponding to tails. $\endgroup$ – David Cary Jul 6 '18 at 18:42
2
$\begingroup$

It appears we are allowed to change the mapping of the outcome of each flip, every time we flip. So, using for convenience the first seven positive integers, we give the following orders:

1st Flip, map $H \to 1$
2nd Flip, map $H \to 2$
...
7th flip, map $H \to 7$
8th flip, map $H \to 1$

etc

Repeat, always in batches of 7 flips. Map the $T$ outcomes to nothing.

SOME REMARKS ON EFFICIENCY
Our RNG, index it by $AP$, will generate zero useful outcomes in one 7-flip batch if we get $T$ in all 7 flips. So

$$P_{AP}(\text{no integers generated}) = (1-p)^7$$

As we run $N_b$ 7-flip batches, the total number of useless flips will tend to

$$\text{Count}_{AP}(\text{useless flips}) \to 7\cdot N_b(1-p)^7$$

Consider now the RNG of @DilipSarwate. There, we use a binomial $B(p,n=5)$ and 5-flip batches. The seven outcomes that generate an integer each has probability of occuring $p^3(1-p)^2$, so, in a 5-flip batch

$$P_{DS}(\text{no integers generated}) = 1-7\cdot p^3(1-p)^2$$

The count of useless flips will here tend to

$$\text{Count}_{DS}(\text{useless flips}) \to 5\cdot N_b\cdot [1-7\cdot p^3(1-p)^2]$$

For the $AP$ RNG to tend to produce less useless flips, it must be the case that

$$ \text{Count}_{AP}(\text{useless flips}) < \text{Count}_{DS}(\text{useless flips})$$

$$\implies 7\cdot N_b(1-p)^7 < 5\cdot N_b\cdot [1-7\cdot p^3(1-p)^2]$$

$$\implies 7\cdot (1-p)^7 < 5\cdot [1-7\cdot p^3(1-p)^2]$$

Numerical examination shows that if $p>0.0467$, then the $AP$ RNG generates less useless flips.

We also find that the number of useless flips is monotonically decreasing in $p$ for the $AP$ RNG, while for the $DS$ RNG it has a minimum at around $p\approx 0.5967$ and then increases again, while in general it stays high. The ratio

$$\frac{\text{Count}_{AP}(\text{useless flips})}{\text{Count}_{DS}(\text{useless flips})}$$

declines pretty quickly. For example it is equal to $0.67$ for$p=0.1$, equal to $0.3$ for $p=0.2$, equal to $0.127$ for $p=0.4$.

$\endgroup$
  • $\begingroup$ Suppose we want to minimize the total number of flips required to generate an integer. I think that your answer provides the minimum total flips across all $p\in(0,1)$. Can you elaborate on that point? $\endgroup$ – Sycorax Jul 6 '18 at 15:58
  • 1
    $\begingroup$ @Sycorax I added something, although I am not sure it is along the lines you suggested. $\endgroup$ – Alecos Papadopoulos Jul 6 '18 at 18:40
  • $\begingroup$ @MartijnWeterings Why? if we get all $H$, we will get each integer once. Why is that a problem? $\endgroup$ – Alecos Papadopoulos Jul 6 '18 at 19:10
  • 1
    $\begingroup$ If you go into this interpretation of (extremely pseudo) "random"-number generator, as just any sequence with equal amount of 1,2,3,4,5,6,7, then you can just as well not flip any coin and directly take the sequence of digits in the number $\frac{1234567}{9999999}$. $\endgroup$ – Martijn Weterings Jul 6 '18 at 20:32
  • 1
    $\begingroup$ This sequence can clearly be distinguished from a sequence of i.i.d uniforms on {1..7}. For all but one specific value of $p$, the probability that the first value will be 1 is not $\frac{1}{7}$. For any value of $p$, the event that the second value is 2 will not be independent of the event that the first value is 1. A sequence which has, on average, an equal number of each value, is a much much weaker condition that i.i.d uniforms. Eg., a sequence where I choose a number fairly and repeat it a million times, satisfies the first condition (and is much more 'efficient' than this answer). $\endgroup$ – jwg Jul 11 '18 at 11:27

protected by kjetil b halvorsen Jul 7 '18 at 12:31

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

Not the answer you're looking for? Browse other questions tagged or ask your own question.