Python for Bayesian Data Analysis

Bayesian data analysis is a statistical paradigm in which uncertainties are modeled as probability distributions rather than single-valued estimates. Because of the ability of its methods to handle uncertain information in a probabilistic fashion, bayesian data analysis methods have become a central concept in data science processes, hence the importance of learning how to implement such methods programatically. The omnipresence of uncertainty in every aspect of our daily lives, including business contexts, also underlines the importance of acquiring Bayesian data analysis skills.

This tutorial illustrates the Python-based application of Bayesian data analysis principles to estimate the average monthly number of tourists visiting the island of Taiwan, based on synthetic data.

Step-by-Step Process

We start by importing the necessary libraries and packages. While there are specialized Python libraries for implementing Bayesian data analysis processes like PyMC3 or cmdstanpy, we keep some degree of simplicity by using the scipy.stats statistical package, still valid for basic bayesian inference processes.

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

In both Bayesian and frequentist statistical methods, data or evidence collection is usually the first part of the process. Therefore, the first instructions will focus on generating some synthetic data, encompassing a sample of 30 instances describing monthly tourist numbers in Taiwan across a two-year and half period (30 months). The data are randomly generated by following a normal distribution with a mean of 10.000 and a standard deviation of 2.000:

np.random.seed(42)
true_mean = 10000
data = np.random.normal(loc=true_mean, scale=2000, size=30)

Below is an example of what a randomly generated sample looks like. Decimals could be optionally removed to reflect better that these values represent people’s numbers.

[10993.42830602  9723.47139766 11295.3770762  13046.05971282
  9531.69325055  9531.7260861  13158.42563101 11534.86945831
  9061.05122813 11085.12008717  9073.16461438  9068.54049286
 10483.92454313  6173.43951068  6550.16433497  8875.42494152
  7974.33775933 10628.49466519  8183.95184896  7175.39259733
 12931.29753784  9548.44739903 10135.05640938  7150.50362757
  8911.23455095 10221.84517942  7698.01284516 10751.39603669
  8798.72262016  9416.61250041]

Next, using parameters -mean and variance- we define a prior probability distribution of the phenomenon to analyze: the mean number of monthly tourists. The initial assumption here (before analyzing the data) is a mean of 0 and a variance of 10000².

A prior distribution represents the a priori belief about the phenomenon studied, which will later be combined with evidence-based likelihoods to yield a posterior distribution.

prior_mean = 0
prior_var = 10000**2
prior_precision = 1 / prior_var

Precision, given by the reciprocal of variance, will also be needed in subsequent stages of the process.

Likelihoods are calculated in parametric form, namely computing the sample mean and variance from the data. Since the generated data are stored in an array, both parameters can be easily obtained using the NumPy library. The ddof argument is set to 1 to calculate an unbiased variance estimate.

data_mean = np.mean(data)
data_var = np.var(data, ddof=1)
n = len(data)
data_precision = 1 / data_var

The above-calculated likelihood parameters explain how likely the observed data values are. The central step in Bayesian data analysis is to combine the likelihood with the prior distribution to obtain an updated distribution, called posterior distribution. This is a fundamental difference concerning frequentist statistical methods, where likelihoods are used to compute p-values or confidence intervals, without combining them with prior distributions.

The following code calculates the posterior mean and posterior variance using conjugate prior. Concretely, prior precision and likelihood precision are summed here to obtain a posterior precision. This is where the initial belief about monthly tourist numbers is updated based on the observed data by forming a posterior distribution.

posterior_precision = prior_precision + n * data_precision
posterior_mean = (prior_precision * prior_mean + n * data_precision * data_mean) / posterior_precision
posterior_var = 1 / posterior_precision

Depending on the random data generated earlier, by printing the posterior estimates we would get something similar to:

Posterior Mean: 9613.32
Posterior Variance: 107885.03

Interestingly, the posterior mean closely approximates the “true mean”, or the mean we initially considered for the data generation.

We finalize by plotting the data sample in an interval histogram along with the posterior mean and also representing the posterior distribution in a separate plot.

plt.figure(figsize=(10, 5))
plt.hist(data, bins=10, alpha=0.75, color='blue', label='Data')
plt.axvline(posterior_mean, color='red', linestyle='dashed', linewidth=1.5, label='Posterior Mean')
plt.title('Histogram of Monthly Tourist Numbers')
plt.xlabel('Number of Tourists')
plt.ylabel('Frequency')
plt.legend()
plt.show()

x = np.linspace(posterior_mean - 3*np.sqrt(posterior_var), posterior_mean + 3*np.sqrt(posterior_var), 100)
posterior_pdf = stats.norm.pdf(x, loc=posterior_mean, scale=np.sqrt(posterior_var))

plt.figure(figsize=(10, 5))
plt.plot(x, posterior_pdf, color='green', label='Posterior Distribution')
plt.title('Posterior Distribution of Tourist Numbers')
plt.xlabel('Number of Tourists')
plt.ylabel('Density')
plt.legend()
plt.show()

Output: