Bayesian Updating in Python

Introduction

In one of my previous posts I discuss how you can update your beliefs using Bayesian updating. You can check out that post here:

Bayesian Updating Simply Explained

In that post, we had three dice with different number ranges. We picked up one die random and did two successive rolls with that die. From this information, we then computed the probability (posterior) of which dice we most likely picked up.

This process was all carried out by hand as we luckily only had two roll outcomes and three dice. However, with more dice and more rolls this problem quickly becomes too long and tedious to do on pen and paper.

In this article, we will implement Bayesian updating for the above problem in Python to expedite the process of calculating the posterior.

Bayes’ Theorem Recap

Lets quickly recap over Bayes’ theorem and its key features:

• P(H): probability of the hypothesis, this is the prior. This is how likely our hypothesis is before we see our data, D.
• P(D|H): the likelihood, the probability of our data being correct given our hypothesis.
• P(H|D): the probability our hypothesis is true from our given data. This is the posterior.
• P(D): the probability of the observed data. This is the normalising constant, that is the sum of the product from the likelihoods and priors:

For a full derivation and intuition of Bayes’ theorem, check out my previous post on the subject:

Conditional Probability and Bayes’ Theorem Simply Explained

The Problem

Lets say have a set of dice with different number ranges, 1–2, 1-3, 1-4 etc. We choose one of these die at random and roll it several times. Using the data we receive from the outcome of the rolls, we can update our belief to estimate which die we most likely picked up.

The Code

Lets start by importing our packages:

import numpy as np
import matplotlib.pyplot as plt

Dice

Generate our dice, where the number represents the highest value on that die. For example, 9 means the die has numbers between 1–9:

dice = np.arange(3,13)
dice

Output: array([ 3,  4,  5,  6,  7,  8,  9, 10, 11, 12])

Prior

As each die is equally likely to be chosen, we have a uniform prior. Therefore, each die has an equal prior probability:

prior = [1/len(dice) for _ in np.arange(len(dice))]
prior

Output: [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]

As we have 10 dice, each die has a probability of 10% of being chosen.

Likelihood and Roll 1

We choose one of the die, and we roll a 1. The likelihoods of rolling a 1 for our dice is as follows:

roll = 1
likelihood1 = [1/i if i >= roll else 0 for i in dice]
likelihood1

Output: [0.3333333333333333,
 0.25,
 0.2,
 0.16666666666666666,
 0.14285714285714285,
 0.125,
 0.1111111111111111,
 0.1,
 0.09090909090909091,
 0.08333333333333333

Notice that the die with 1–3 has the highest likelihood. This makes sense as it has the smallest number range.

Computing Posterior

Multiplying the priors and likelihoods, we can find the posterior and normalise it:

posterior = np.array(likelihood1) * np.array(prior)
list(posterior/sum(posterior))

Ouput: [0.20791611349879613,
 0.15593708512409712,
 0.1247496680992777,
 0.10395805674939806,
 0.08910690578519834,
 0.07796854256204856,
 0.06930537116626538,
 0.06237483404963885,
 0.05670439459058077,
 0.05197902837469903]

Plotting the posterior:

plt.figure(figsize=(13,7))
plt.xlabel('Dice', fontsize=20)
plt.ylabel('Probability', fontsize=20)
plt.plot(dice, list(posterior/sum(posterior)))
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.show()

Finding the die with the highest probability, this is known as the maximum a posteriori probability (MAP):

dice[np.argmax(posterior)]

Output: 3

Therefore, the most likely die is the one with ranges from 1–3! This is quite obvious since it had the highest likelihood and we also had a uniform prior.

Likelihood and Roll 2

We roll the same die again and get a 5. The likelihood of this outcome is:

roll = 5
likelihood2 = [1/i if i >= roll else 0 for i in dice]
likelihood2

Output: [0,
 0,
 0.2,
 0.16666666666666666,
 0.14285714285714285,
 0.125,
 0.1111111111111111,
 0.1,
 0.09090909090909091,
 0.08333333333333333]

Notice the first two are now 0. This is because it is impossible to roll a 5 with dice that have ranges 1–3 and 1–4.

Update Posterior

Using our old posterior as the new prior along with the likelihoods of rolling a 5, we update our posterior as:

posterior = posterior * np.array(likelihood2)
list(posterior/sum(posterior))

Output: [0.0,
 0.0,
 0.2829544144262495,
 0.1964961211293399,
 0.14436449715624972,
 0.1105290681352537,
 0.08733160939081774,
 0.07073860360656238,
 0.05846165587319205,
 0.049124030282334974]

Plotting the posterior:

plt.figure(figsize=(13,7))
plt.xlabel('Dice', fontsize=20)
plt.ylabel('Probability', fontsize=20)
plt.plot(dice, list(posterior/sum(posterior)))
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.show()

The die with the highest probability:

dice[np.argmax(posterior)]

Output: 5

The die with the range 1–5 is the most likely given our data!

Dice with ranges 1–3 and 1–4 are 0 because it is impossible for them to output a 5, hence we know we definitely didn’t pick them up!

Generic Function

Lets wrap all this code into a function for convenience:

def bayesian_dice_updating(data, dice):
    """ Compute the posterior distribution for given dice and data.

    :param data: The numbers that have been rolled from the dice
    :type data: list, np.array

    :param dice: The range of dices where the number represents 
                 the maximum value that die can take.
    :type dice: list, np.array

    :returns: Posterior distribution of the dice given the data
    :rtype: list

    """

    prior = [1/len(dice) for _ in np.arange(len(dice))]

    posterior = prior

    for roll in data:
        likelihood = [1/i if i >= roll else 0 for i in dice]
        posterior = np.array(likelihood) * np.array(posterior)

    return list(posterior/sum(posterior))

Let’s try it out!

Generate some dice and some data (roll outcomes):

dice = np.arange(1,51)
data = [4,6,9]

Plot the results:

plt.figure(figsize=(13,7))
plt.xlabel('Dice', fontsize=20)
plt.ylabel('Probability', fontsize=20)
plt.plot(dice, bayesian_dice_updating(data, dice))
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.show()

dice[np.argmax(bayesian_dice_updating(data, dice))]

Output: 9

The functions works well!

The full code is available at my GitHub here:

Medium-Articles/Bayesian_Updating.ipynb at main · egorhowell/Medium-Articles

Conclusion

In this article we recapped over Bayes’ theorem and showed how to code up Bayesian updating in Python to make computing the posterior easy for a simple example.

Another Thing!

I have a free newsletter, Dishing the Data, where I share weekly tips for becoming a better Data Scientist. There is no "fluff" or "clickbait," just pure actionable insights from a practicing Data Scientist.

Dishing The Data | Egor Howell | Substack

Connect With Me!

• YouTube
• LinkedIn
• Twitter
• GitHub