Bayesian updating with PYMC

I’ve been looking for neat ways to update a Bayesian Prior from a posterior sample for a while, and just the other day managed to find what I was looking for: a code example that shows how to make a posterior to a distribution, known to PYMC. The code for doing that, by jcrudy, can be found here. The key insight with that code is the “known to PYMC”-part: it’s no problem to make a posterior sample to a distribution, but in order to make it compatible with PYMC, it needs to be done as a PYMC object.  And that’s what jcrudy’s code does.

Anyway, below a small example, borrowed from Richard McElreath’s excellent book “Statistical Rethinking”.

Let’s pretend that we wanted to know the ratio of Water to Land of our earth, and that the way we will go about it is to toss an inflatable globe into air a number of times, and keep track of where we first touch it when we catch it, land or water.

Let’s say that after the first 9 tosses, our results look like below:

[1,0,1,0,1,0,1,1,1] where ‘1’ indicates water, and ‘0’ land.

The basic logic of our Bayesian Inference is to update the prior,  each time we get a new data point, i.e. either a ‘1’ or a ‘0’.

Let’s furthermore say that we have an initial idea of the correct ratio being somewhere around 50%, but we are far from certain, so we will set a fairly “loose” initial prior, a Beta-distribution centered on 50%, as below. This prior allows the calculated rate to be almost the entire range of 0% to 100%, with the highest probability around 50%.

Now, what we want to do is to update our prior belief – the one initially being centered around 50% – based on each subsequent data point.

Below the results of 18 of such updates (the above vector of 9 points is concatenated to provide these 18 data points):

This graph is read from left-to-right, and shows the updated prior, or current belief on the correct ration, after each subsequent data point has been processed.

From the topmost left image, we can see that when the first data point, which happens to be a ‘1’, has been processed, the Bayesian model has adjusted its belief on the correct ratio somewhat to the right, from 50% to 55%, and after having processed the ‘0’ in the second data point, reduced the current belief back to 50%. And so on so forth as the remaining data points get processed.

Finally, after having processed all 18 data points, updating the prior for each point, the Bayesian model settles for a ratio of 64% (mean value). But the real benefit of Bayesian analysis is that we actually get a probability distribution as a result, not just a point estimate. Thus, the Bayesian approach preserves the underlying uncertainty of the model in a very neat, visually obvious way.

If you look closely at the 18 different plots above, you will see that the distribution gets narrower as each data point is processed. This means the Bayesian model is gradually getting more and more certain about the sought after ratio, as it gets more data to base its estimate on.

Code below.

import prior_kde as pk
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as sps

data = [1,0,1,0,1,0,1,1,1] * 2

# Initial Prior
prior = pm.Beta('prior',alpha=5,beta=5)

i = 0
cols = 3
rows = 6

plt.figure(figsize=(18,12))
plt.title('Initial Prior')
initial_prior_samples = [prior.random() for i in range(10000)]
plt.hist(initial_prior_samples)
plt.axvline(np.mean(initial_prior_samples),color='r',ls='dashed')
plt.savefig('initial_prior.jpg',format='jpg')
           
fig,axes = plt.subplots(rows,cols,sharex=True,sharey=True,
                        figsize=(18,12))

for d in data:
    print ('\n***** iteration {} *****'.format(i+1))

    # Likelihood
    obs = pm.Bernoulli('obs',p=prior,value=d,observed=True)

    model = pm.Model([prior,obs])
    mcmc = pm.MCMC(model)
    mcmc.sample(40000,10000)

    prior_samples = np.random.choice(mcmc.trace('prior')[:],replace=True,
                                     size=5000)

    mean = np.mean(prior_samples)
    CI = np.percentile(prior_samples,[2.5,97.5])
    
    row = i // cols  
    col = i % cols 

    bars = axes[row,col].hist(prior_samples,bins=30)
    
    axes[row,col].set_title(str('{},  mean: {:.2}').format(
        str(data[0:i+1]),np.mean(prior_samples)))
    axes[row,col].axvline(np.mean(prior_samples),color='r',ls='dashed')
    
    # Update Prior by making a pymc object of samples
    prior = pk.KernelSmoothing('prior',prior_samples)

    i += 1

plt.tight_layout()
plt.savefig('updated_priors.jpg',format='jpg')
plt.show()