# Rule of Succession for Unfair Coin

Given the first n flip results from an unfair coin, we wish to estimate the probability that the next flip is a heads. I can take 2 approaches to this:

Frequentist: P(Next is heads|previous sequence of flips) = p (independence), and the MVUE for p is p_hat_freq = #heads/#flips

Bayesian: Assume a uniform(0,1) prior for p, then condition the calculate the posterior P(Next is heads|previous sequence, H), where H is our prior hypothesis, ..., giving the well-known Laplace Rule of Succession p_hat_Bayes = (#heads + 1)/(#flips + 2). For a complete derivation see https://en.wikipedia.org/wiki/Rule_of_succession

If I place myself at that decision point, I can use either of the two estimators above. When simulating this situation, the frequentist approach always wins, even when using a uniform prior for selecting p. Let's say we draw p, flip 10 times, and use the two estimators to estimate the true p. The frequentist estimator is unbiased, while the Bayesian is only asymptotically unbiased, and is biased for finite samples. Here is a sample run, along with the code.

My question is - faced with the decision after seeing the first 10 flips, as in the simulation, why would I ever use the Bayesian estimate? It is biased, performs worse empirically, and even has lower standard error. So it is biased, but more sure of itself. The only benefit I see is when the number of flips and number of trials are both very low, in the ranges of 1-10 for each.

from __future__ import division
import numpy as np
import numpy.random as rnd
import pandas as pd
import matplotlib.pyplot as plt
import pdb

def calc_probs(p,num_flips,num_trials):

freq_est  = []
Bayes_est = []

for trial_num in np.arange(num_trials):
devs = rnd.uniform(0,1,num_flips)
flips = [1*(x<p) for x in devs]

# Frequentist estimate

# Bayesian estimate

return pd.DataFrame({ 'num_trials' : num_trials,
'p'          : [p],
'Bayes_est'  : [np.mean(Bayes_est)],
'freq_est'   : [np.mean(freq_est)],
'Bayes_sd'   : [np.sqrt(np.var(Bayes_est))],
'freq_sd'    : [np.sqrt(np.var(freq_est))]})

if __name__ == "__main__":

if (True):
num_flips = 10
num_trials = 100
r = pd.DataFrame()
for num_trial in np.arange(10):
p = rnd.uniform(0,1)
r_loc = calc_probs(p,num_flips,num_trials)
r_loc['trial_num'] = 1+num_trial
r = pd.concat([r,r_loc])

fig,ax = plt.subplots(nrows=1,ncols=1)
ax.errorbar(r['trial_num'],r['Bayes_est'],yerr=0,fmt='b*-',label='Bayes')
ax.errorbar(r['trial_num'],r['freq_est'],yerr=0,fmt='r*-',label='Freq')
ax.plot(r['trial_num'], r['p'],'k--',linewidth=2,label='Actual')
ax.set_title('Loaded Coin, Num_flips: ' + str(num_flips))
ax.set_xlabel('Trial Number')
ax.set_ylabel('p_hat')
ax.legend(loc='upper right')
plt.show()

• There is neither one Bayesian nor one Classical ("Frequentist") estimator. There are infinitely many of both--and every Bayesian point estimator of the probability is also a Classical estimator. However, there can be Classical estimators that are non-Bayes; and Bayesian estimators (of course) consist of the entire posterior distribution. These considerations suggest your dichotomy does not correctly capture any inherent differences between Bayesian and Classical procedures.
– whuber
Commented May 22, 2016 at 22:29

I think you are missing the case of events where you don't observe a success/failure within the first 10 flips. Your Frequentist estimation would be that the probability of say a Heads would be either 1 or 0, which you may know to be untrue. The Bayesian approach would leave the possibility as finite though unlikely. This is discussed in the Rule of Succession link you provided.

If you want to argue that there are such possibilities where it is impossible to get a success or failure then I would suggest that you are using an improper prior. I would suggest you use the prior discussed in your linked Rule of Succession article that gives a better prior for that case.

I should also note that your metric for success is kind of weird in my opinion. The situation you give is having just flipped 10 coins and estimating the probability on those flips. I don't know how you you want to treat being off, whether it is better to underestimate, overestimate, or if both equally bad, but your metric currently doesn't do that. It instead takes the average of the predictions and packages them, this would violate the one universe idea.

I took your code and ported it to R since I'm more comfortable with that language and it has some built in features that are perfect for what we are trying to do. I also added a mean absolute error metric (quick and easy to implement):

set.seed(10)

calc_probs <- function(num_flips,num_trials)
{
p <- runif(1,0,1)
results <- rbinom(n = num_trials,size = num_flips, prob = p)
freq_est <- sapply(results, function(x) x/num_flips)
bayes_est <- sapply(results, function(x) (x+1)/(num_flips+2))
return(data.frame(p = p,
freq_est = mean(freq_est),
bayes_est = mean(bayes_est),
freq_var = sqrt(var(freq_est)),
bayes_var = sqrt(var(bayes_est)),
freq_error = mean(abs(freq_est-p)),
bayes_error = mean(abs(bayes_est-p))
)
)
}

num_flips <- 10
num_trials <- 100
num_runs <- 10000
scale <- 10

temp <- calc_probs(num_flips,num_trials)
for(i in 1:(num_runs-1))
temp <- rbind(temp, calc_probs(num_flips,num_trials))

temp$binned_p <- floor(temp$p*scale)/scale
temp$bayes_better <- temp$bayes_error<temp$freq_error bayes_better_quantiles <- aggregate(bayes_better~binned_p,temp,mean) freq_error_quantiles <- aggregate(freq_error~binned_p,temp,mean) bayes_error_quantiles <- aggregate(bayes_error~binned_p,temp,mean) hist(temp$p)
plot(freq_error_quantiles)
plot(bayes_error_quantiles)
plot(bayes_better_quantiles)

print(paste0("Average frequentist estimate error: ", mean(temp$freq_error))) print(paste0("Average bayesian estimate error: ", mean(temp$bayes_error)))


Starting with a histogram of $p$: Next we examine the error by decile for the frequentist estimate: Finally we look at the error by decile for the bayesian estimate you provided: Comparing for when the bayesian approach is better (note this isn't by how much so it is a bit of a cheat): We notice that the frequentist estimate has less average error when $p$ is more extreme (near 0 or 1) and the bayes approach when $p$ is closer to 1/2. This is expected based on the prior we chose for the bayes approach. Overall the average error calculated by my code is:

"Average frequentist estimate error: 0.0999418570390359"
"Average bayesian estimate error: 0.0928987467885039"


So the average mean absolute error is better for the Bayesian approach.

One may point out that I originally claimed that the point of the Rule of Succession was to avoid problems when the probability was low or high. In order to better demonstrate the ability of the Bayesian Method I would like to use a more standard error metric for these types of problems. In particular the log loss function. This can be done by adjusting my calc_probs function as follows (If this is wrong someone point this out so I can fix it.):

calc_probs <- function(num_flips,num_trials, cutoff = 0.01)
{
p <- runif(1,0,1)
results <- rbinom(n = num_trials,size = num_flips, prob = p)
freq_est <- sapply(results, function(x) x/num_flips)
bayes_est <- sapply(results, function(x) (x+1)/(num_flips+2))
return(data.frame(p = p,
freq_est = mean(freq_est),
bayes_est = mean(bayes_est),
freq_var = sqrt(var(freq_est)),
bayes_var = sqrt(var(bayes_est)),
freq_error = mean(abs(freq_est-p)),
bayes_error = mean(abs(bayes_est-p)),
freq_log_odds_error = mean(p*(log(max(freq_est,cutoff)))+
(1-p)*log(max(cutoff,1-freq_est))),
bayes_log_odds_error = mean(p*log(max(bayes_est,cutoff))+
(1-p)*log(max(cutoff,1-bayes_est)))
)
)
}


Doing so running my code to plot these loss functions I get the following results.

For the frequentist method the binned mean log-loss function is:

For the bayesian method I get:

For the which is better comparison I get:

Finally overall loss I get:

"Average frequentist estimate log loss: -0.110217675806207"
"Average bayesian estimate log loss -0.184244915921267"


A reasonable objection to my implementation of log-loss here is that I have used $p$ instead of making an eleventh prediction then doing a classification error based on that. This is a reasonable objection and I will have to return to update my method to take this into account.

Ignoring the objection we see that universally the Bayesian method works better for the log-loss error function. This is inline with the Rule of Succession idea since its purpose is to avoid a $p=0$ prediction for events that are possible.

My question is - faced with the decision after seeing the first 10 flips, as in the simulation, why would I ever use the Bayesian estimate? It is biased, performs worse empirically, and even has lower standard error.

I think the problem here is that the simulation isn't properly capturing the situation in which you have to make a decision after seeing ten flips.

The question that both estimators are trying to answer is: Given that I flipped a coin ten times and got a certain number of heads, what is my best guess for $p$?

What you did is: simulate this process (of flipping a coin ten times and getting an estimate for $p$) one hundred times, and then averaging the estimates for $p$ over these hundred repeats. But the average of these estimates isn't in general an estimate for $p$, unless you are treating the whole process as flipping a single coin $1000$ times, in which case your Bayesian estimate would be ((number of heads) +1)/1002.

Since the frequentist estimate of $p$ is unbiased, and you used the same value of $p$ for every trial, it's absolutely to be expected that the average of $100$ of these estimates will be close to the true $p$, because this is the definition of unbiasedness. Since the Bayesian estimate is biased, in general averaging $100$ of them is not going to get you closer to the true $p$.

What you should do to compare the two estimators of $p$ is: repeatedly simulate ten coin flips, compute the two estimates of $p$ and see which one is closer to the true $p$. In other words, you should look at something like

$$\text{average over many trials of } | \text{true } p - \text{estimated } p |$$

instead of what you did, which was

$$| \text{true } p - \text{average over many trials of estimated } p |.$$

You will find that the Bayesian estimate is better, unless you pick the $p$'s to be close to $0$ or $1$ most of the time (for example, if you take $p=0.9$ for every trial, then the frequentist estimate will be better.)

• Yes - it looks like you're correct. I have calculated your second expression. Tied up now but I will test tonight. Looks like I can just set num_trials = 1, then do some averaging outside the function. Commented May 23, 2016 at 21:29