# Revisiting the coin-flipping problem

Hello,

Imagine we tossed a biased coin 8 times (we don’t know how biased it is), and we recorded 5 heads (H) to 3 tails (T) so far. What is the probability of that the next 3 tosses will all be tails? In other words, we are wondering the expected probability of having 5Hs and 6Ts after 11th tosses.

I want to build a MCMC simulation model using pyMC3 to find the Bayesian solution. There is also an analytical solution within the Bayesian approach for this problem. So, I will be able to compare the results derived from the simulation, the analytical way as well as the classical frequentest way. Let me briefly explain what I could do so far:

1. Frequentist solution:
If we consider the probability for a single toss: E(T) = p = (3/8) = 0,375
Then, the ultimate answer is E({T,T,T}) = p^3 = (3/8)^3 = 0,052.

2.1. Bayesian solution with analytical way:
Please assume the unknown parameter “p” represents for the bias of the coin.
If we consider the probability for a single toss:
E(T) = Integral0-1[p * P(p | H = 5, T = 3) dp] = 0,400 (I calculated the result after some algebraic manipulation)
E({T,T,T}) = Integral0-1[p^3 * P(p | H = 5, T = 3) dp] = 10/11 = 0,909.

2.2. Bayesian solution with MCMC simulation:
When we consider the probability for a single toss, I built the model in pyMC3 as below:

Tail: 1

``````data = [0, 0, 0, 0, 0, 1, 1, 1]
import pymc3 as pm
with pm.Model() as coin_flipping:
p = pm.Uniform(‘p’, lower=0, upper=1)
y = pm.Bernoulli(‘y’, p=p, observed=data)
trace = pm.sample(1000)
pm.traceplot(trace)
``````

After I run this code, I got that the posterior mean is E(T) =0,398 which is very close to the result of analytical solution (0,400). I am happy so far, but this is not the ultimate answer. I need a model that simulate the probability of E({T,T,T}). I appreciate if someone help me.

You can do it like this:

``````import numpy as np
import theano
import pymc3 as pm
data = theano.shared(np.asarray([0, 0, 0, 0, 0, 1, 1, 1]))

with pm.Model() as coin_flipping:
p = pm.Uniform('p', lower=0, upper=1)
y = pm.Bernoulli('y', p=p, observed=data)
trace = pm.sample(1000, tune=1000)
# reset value to get the shape right
data.set_value(np.asarray([0, 0, 0]))
ppc = pm.sample_posterior_predictive(trace)

np.mean(ppc['y'].sum(axis=1)==3)
``````
1 Like

Many thanks Junpeng,
But, It raised that: AttributeError: module ‘pymc3’ has no attribute ‘sample_posterior_predictive’
What should I do?

1 Like

Thanks Junpeng, that’s really wonderful

@junpenglao in your example, why did you do this?

``````# reset value to get the shape right
data.set_value(np.asarray([0, 0, 0]))
``````

I thought 0 indicated a Head while 1 indicated a Tail. @serdar asked “What is the probability of that the next 3 tosses will all be tails?”, so shouldn’t you set the value to [1,1,1]?

`sample_posterior_predictive` will convert a trace (in this case over `p`) to the posterior predictive distribution over `y`.

This distribution is `k` independent draws of a bernoulli, with probability `p`; and `k=data.shape[0]`.

Notably, because we are by design sampling over the entire probability space, only the shape of data matters, and not its value. You could even use

`data.set_value(np.as_array([-4,-4,-4]))`

and shouldn’t (in theory) see any complaint. Computation of `p[TTT]` happens after.