Revisiting the coin-flipping problem

Questions
1

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)
Similarly, the ultimate answer is:
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:

Head: 0
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.

2

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
3

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

4

Please try upgrading your PyMC3 :wink:

1 Like
5

Thanks Junpeng, that’s really wonderful :blush:

6

@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]?

7

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.