### Prerequisites

This post is an introduction to Bayesian probability and inference. We will discuss the intuition behind these concepts, and provide some examples written in Python to help you get started. To get the most out of this introduction, the reader should have a basic understanding of statistics and probability, as well as some experience with Python. The examples use the Python package pymc3.

### Introduction to Bayesian Thinking

Bayesian inference is an extremely powerful set of tools for modeling any random variable, such as the value of a regression parameter, a demographic statistic, a business KPI, or the part of speech of a word. We provide our understanding of a problem and some data, and in return get a quantitative measure of how certain we are of a particular fact. This approach to modeling uncertainty is particularly useful when:

• Data is limited
• We have reason to believe that some facts are more likely than others, but that information is not contained in the data we model on
• We're interested in precisely knowing how likely certain facts are, as opposed to just picking the most likely fact

The table below enumerates some applied tasks that exhibit these challenges, and describes how Bayesian inference can be used to solve them. Don't worry if the Bayesian solutions are foreign to you, they will make more sense as you read this post:

Typically, Bayesian inference is a term used as a counterpart to frequentist inference. This can be confusing, as the lines drawn between the two approaches are blurry. The true Bayesian and frequentist distinction is that of philosophical differences between how people interpret what probability is. We'll focus on Bayesian concepts that are foreign to traditional frequentist approaches and are actually used in applied work, specifically the prior and posterior distributions.

Consider Bayes' theorem:

Think of A as some proposition about the world, and B as some data or evidence. For example, A represents the proposition that it rained today, and B represents the evidence that the sidewalk outside is wet:

p(rain | wet) asks, "What is the probability that it rained given that it is wet outside?" To evaluate this question, let's walk through the right side of the equation. Before looking at the ground, what is the probability that it rained, p(rain)? Think of this as the plausibility of an assumption about the world. We then ask how likely the observation that it is wet outside is under that assumption, p(wet | rain)? This procedure effectively updates our initial beliefs about a proposition with some observation, yielding a final measure of the plausibility of rain, given the evidence.

This procedure is the basis for Bayesian inference, where our initial beliefs are represented by the prior distribution p(rain), and our final beliefs are represented by the posterior distribution p(rain | wet). The denominator simply asks, "What is the total plausibility of the evidence?", whereby we have to consider all assumptions to ensure that the posterior is a proper probability distribution.

Bayesians are uncertain about what is true (the value of a KPI, a regression coefficient, etc.), and use data as evidence that certain facts are more likely than others. Prior distributions reflect our beliefs before seeing any data, and posterior distributions reflect our beliefs after we have considered all the evidence. To unpack what that means and how to leverage these concepts for actual analysis, let's consider the example of evaluating new marketing campaigns.

### Example: Evaluating New Marketing Campaigns Using Bayesian Inference

Assume that we run an ecommerce platform for clothing and in order to bring people to our site, we deploy several digital marketing campaigns. These campaigns feature various ad images and captions, and are presented on a number of social networking websites. We want to present the ads that are the most successful. For the sake of simplicity, we can assume that the most successful campaign is the one that results in the highest click-through rate: the ads that are most likely to be clicked if shown.

We introduce a new campaign called "facebook-yellow-dress," a campaign presented to Facebook users featuring a yellow dress. The ad has been presented to 10 users so far, and 7 of the users have clicked on it. We would like to estimate the probability that the next user will click on the ad.

By encoding a click as a success and a non-click as a failure, we're estimating the probability θ that a given user will click on the ad. Naturally, we are going to use the campaign's historical record as evidence. Because we are considering unordered draws of an event that can be either 0 or 1, we can infer the probability θ by considering the campaign's history as a sample from a binomial distribution, with probability of success θ. Traditional approaches of inference consider multiple values of θ and pick the value that is most aligned with the data. This is known as maximum likelihood, because we're evaluating how likely our data is under various assumptions and choosing the best assumption as true. More formally:

argmaxθp(X |θ), where X is the data we've observed.

Here, p(X |θ) is our likelihood function; if we fix the parameter θ, what is the probability of observing the data we've seen? Let's look at the likelihood of various values of θ given the data we have for facebook-yellow-dress:

import numpy as np
from scipy.misc import factorial
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (16,7)

def likelihood(theta, n, x):
"""
likelihood function for a binomial distribution

n: [int] the number of experiments
x: [int] the number of successes
theta: [float] the proposed probability of success
"""
return (factorial(n) / (factorial(x) * factorial(n - x))) \
* (theta ** x) * ((1 - theta) ** (n - x))

#the number of impressions for our facebook-yellow-dress campaignn_impressions = 10.

#the number of clicks for our facebook-yellow-dress campaign
n_clicks = 7.
#observed click through rate
ctr = n_clicks / n_impressions
#0 to 1, all possible click through rates
possible_theta_values = map(lambda x: x/100., range(100))

#evaluate the likelihood function for possible click through rates
likelihoods = map(lambda theta: likelihood(theta, n, x)\
, possible_theta_values)

#pick the best theta
mle = possible_theta_values[np.argmax(likelihoods)]
#plot
f, ax = plt.subplots(1)
ax.plot(possible_theta_values, likelihoods)
ax.axvline(mle, linestyle = "--")
ax.set_xlabel("Theta")
ax.set_ylabel("Likelihood")
ax.grid()
ax.set_title("Likelihood of Theta for New Campaign")
plt.show()

Of the 10 people we showed the new ad to, 7 of them clicked on it. So naturally, our likelihood function is telling us that the most likely value of theta is 0.7. However, some of our analysts are skeptical. The performance of this campaign seems extremely high given how our other campaigns have done historically. Let's overlay this likelihood function with the distribution of click-through rates from our previous 100 campaigns:

plt.rcParams['figure.figsize'] = (16, 7)
import numpy as np
import pandas as pd

true_a = 11.5
true_b = 48.5
#number of marketing campaigns
N = 100#randomly generate "true" click through rate for each campaign
p = np.random.beta(true_a,true_b, size=N)
#randomly pick the number of impressions for each campaign
impressions = np.random.randint(1, 10000, size=N)
#sample number of clicks for each campaign
clicks = np.random.binomial(impressions, p).astype(float)
click_through_rates = clicks / impressions
#plot the histogram of previous click through rates with the evidence#of the new campaign
f, ax = plt.subplots(1)
ax.axvline(mle, linestyle = "--")
ax.plot(possible_theta_values, likelihoods)

zero_to_one = [j/100. for j in xrange(100)]
counts, bins = np.histogram(click_through_rates
, bins=zero_to_one)
counts = counts / 100.
ax.plot(bins[:-1],counts, alpha = .5)
line1, line2, line3 = ax.lines
ax.legend((line2, line3), ('Likelihood of Theta for New Campaign'
, 'Frequency of Theta Historically')
, loc = 'upper left')
ax.set_xlabel("Theta")
ax.grid()
ax.set_title("Evidence vs Historical Click Through Rates")
plt.show()

Clearly, the maximum likelihood method is giving us a value that is outside what we would normally see. Perhaps our analysts are right to be skeptical; as the campaign continues to run, its click-through rate could decrease. Alternatively, this campaign could be truly outperforming all previous campaigns. We can't be sure. Ideally, we would rely on other campaigns' history if we had no data from our new campaign. And as we got more and more data, we would allow the new campaign data to speak for itself.

## The Prior Distribution

This skepticism corresponds to prior probability in Bayesian inference. Before considering any data at all, we believe that certain values of θ are more likely than others, given what we know about marketing campaigns. We believe, for instance, that p(θ = 0.2)>p(θ = 0.5), since none of our previous campaigns have had click-through rates remotely close to 0.5. We express our prior beliefs of θ with p(θ). Using historical campaigns to assess p(θ) is our choice as a researcher. Generally, prior distributions can be chosen with many goals in mind:

• Informative; empirical: We have some data from related experiments and choose to leverage that data to inform our prior beliefs. Our prior beliefs will impact our final assessment.
• Informative; non-empirical: We have some inherent reason to prefer certain values over others. For instance, if we want to regularize a regression to prevent overfitting, we might set the prior distribution of our coefficients to have decreasing probability as we move away from 0. Our prior beliefs will impact our final assessment.
• Informative; domain-knowledge: Though we do not have supporting data, we know as domain experts that certain facts are more true than others. Our prior beliefs will impact our final assessment.
• Non-informative: Our prior beliefs will have little to no effect on our final assessment. We want the data to speak for itself.

For our example, because we have related data and limited data on the new campaign, we will use an informative, empirical prior. We will choose a beta distribution for our prior for θ. The beta distribution is a 2 parameter (α, β) distribution that is often used as a prior for the θ parameter of the binomial distribution. Because we want to use our previous campaigns as the basis for our prior beliefs, we will determine α and β by fitting a beta distribution to our historical click-through rates. Below, we fit the beta distribution and compare the estimated prior distribution with previous click-through rates to ensure the two are properly aligned:

from scipy.stats import beta
#fit beta to previous CTRs
prior_parameters = beta.fit(click_through_rates
, floc = 0
, fscale = 1)#extract a,b from fit
prior_a, prior_b = prior_parameters[0:2]

#define prior distribution sample from prior
prior_distribution = beta(prior_a, prior_b)
#get histogram of samples
prior_samples = prior_distribution.rvs(10000)
#get histogram of samples
fit_counts, bins = np.histogram(prior_samples
, zero_to_one)#normalize histogram
fit_counts = map(lambda x: float(x)/fit_counts.sum()
, fit_counts)
#plot
f, ax = plt.subplots(1)
ax.plot(bins[:-1], fit_counts)

hist_ctr, bins = np.histogram(click_through_rates
, zero_to_one)
hist_ctr = map(lambda x: float(x)/hist_ctr.sum()
, hist_ctr)
ax.plot(bins[:-1], hist_ctr)
estimated_prior, previous_click_through_rates = ax.lines
ax.legend((estimated_prior, previous_click_through_rates)
,('Estimated Prior'
, 'Previous Click Through Rates'))
ax.grid()
ax.set_title("Comparing Empirical Prior with Previous Click Through Rates")
plt.show()

We find that the best values of α and β are 11.5 and 48.5, respectively. The beta distribution with these parameters does a good job capturing the click-through rates from our previous campaigns, so we will use it as our prior. We will now update our prior beliefs with the data from the facebook-yellow-dress campaign to form our posterior distribution.

### The Posterior Distribution

After considering the 10 impressions of data we have for the facebook-yellow-dress campaign, the posterior distribution of θ gives us plausibility of any click-through rate from 0 to 1.

The effect of our data, or our evidence, is provided by the likelihood function, p(X|θ). What we are ultimately interested in is the plausibility of all proposed values of θ given our data or our posterior distribution p(θ|X). From the earlier section introducing Bayes' Theorem, our posterior distribution is given by the product of our likelihood function and our prior distribution:

Since p(X) is a constant, as it does not depend on θ, we can think of the posterior distribution as:

We'll now demonstrate how to estimate p(θ|X) using PyMC.

### Doing Bayesian Inference with PyMC

Usually, the true posterior must be approximated with numerical methods. To see why, let's return to the definition of the posterior distribution:

The denominator p(X) is the total probability of observing our data under all possible values of θ. A more descriptive representation of this quantity is given by:

Which sums the probability of X over all values of θ. This integral usually does not have a closed-form solution, so we need an approximation. One method of approximating our posterior is by using Markov Chain Monte Carlo (MCMC), which generates samples in a way that mimics the unknown distribution. We begin at a particular value, and "propose" another value as a sample according to a stochastic process. We may reject the sample if the proposed value seems unlikely and propose another. If we accept the proposal, we move to the new value and propose another.

PyMC is a python package for building arbitrary probability models and obtaining samples from the posterior distributions of unknown variables given the model. In our example, we'll use MCMC to obtain the samples.

The prototypical PyMC program has two components:

• Define all variables, and how variables depend on each other
• Run an algorithm to simulate a posterior distribution

Let's now obtain samples from the posterior. We select our prior as a Beta(11.5,48.5). Let's see how observing 7 clicks from 10 impressions updates our beliefs:

import pymc3 as pm
import numpy as np

#create our data:clicks = np.array([n_clicks])
#clicks represents our successes. We observed 7 clicks.impressions = np.array([n_impressions])
#this represents the number of trials. There were 10 impressions.

with pm.Model() as model:
#sets a context; all code in block "belongs" to the model object

theta_prior = pm.Beta('prior', 11.5, 48.5)
#our prior distribution, Beta (11.5, 48.5)
observations = pm.Binomial('obs',n = impressions
, p = theta_prior
, observed = clicks)     #Sampling distribition of outcomes in the dataset.
#our prior p_prior will be updated with data

start = pm.find_MAP()    #find good starting values for the sampling algorithm
#Max Aposterior values, or values that are most likely

step = pm.NUTS(state=start)     #Choose a particular MCMC algorithm     #we'll choose NUTS, the No U-Turn Sampler (Hamiltonian)

trace = pm.sample(5000
, step
, start=start
, progressbar=True)               #obtain samples 

Let's walk through each line of code:

with pm.Model() as model:

pm.Model creates a PyMC model object. as model assigns it to the variable name "model", and the with ... : syntax establishes a context manager. All PyMC objects created within the context manager are added to the model object.

theta_prior = pm.Beta('prior', 11.5, 48.5)

Theta_prior represents a random variable for click-through rates. It will serve as our prior distribution for the parameter θ, the click-through rate of our facebook-yellow-dress campaign. This random variable is generated from a beta distribution (pm.Beta); we name this random variable "prior" and hardcode parameter values 11.5 and 48.5. We could have set the values of these parameters as random variables as well, but we hardcode them here as they are known.

observations = pm.Binomial('obs',n = impressions , p = theta_prior , observed = clicks)

This statement represents the likelihood of the data under the model. Again we define the variable name and set parameter values with n and p. Note that for this variable, the parameter p is assigned to a random variable, indicating that we are trying to model that variable. Lastly, we provide observed instances of the variable (i.e. our data) with the observed keyword. Because we have said this variable is observed, the model will not try to change its values.

start = pm.find_MAP()

step = pm.NUTS(state=start)

trace = pm.sample(2000, step, start=start, progressbar=True)

These three lines define how we are going to sample values from the posterior. pm.find_MAP() will identify values of theta that are likely in the posterior, and will serve as the starting values for our sampler.

pm.NUTS(state=start) will determine which sampler to use. The sampling algorithm defines how we propose new samples given our current state. The proposals can be done completely randomly, in which case we'll reject samples a lot, or we can propose samples more intelligently. NUTS (short for the No-U-Turn sample) is an intelligent sampling algorithm. Other choices include Metropolis Hastings, Gibbs, and Slice sampling.

Lastly, pm.sample(2000, step, start=start, progressbar=True) will generate samples for us using the sampling algorithm and starting values defined above.

Let's take the histogram of the samples obtained from PyMC to see what the most probable values of θ are, compared with our prior distribution and the evidence (likelihood of our data for each value of θ):

#plot the histogram of click through rates
plt.rcParams['figure.figsize'] = (16, 7)
#get histogram of samples from posterior distribution of CTRs
posterior_counts, posterior_bins = np.histogram(trace['prior']
,bins=zero_to_one)
#normalized histogramposterior_counts = posterior_counts / float(posterior_counts.sum())
#take the mean of the samples as most plausible value
most_plausible_theta = np.mean(trace['prior'])
#histogram of samples from prior distribution
prior_counts, bins = np.histogram(prior_samples
, zero_to_one)#normalize
prior_counts = map(lambda x: float(x)/prior_counts.sum()
, prior_counts)
#plot
f, ax = plt.subplots(1)
ax.plot(possible_theta_values, likelihoods)
ax.plot(bins[:-1],prior_counts, alpha = .2)
ax.plot(bins[:-1],posterior_counts)
ax.axvline(most_plausible_theta, linestyle = "--", alpha = .2)
line1, line2, line3, line4 = ax.lines
ax.legend((line1, line2, line3, line4), ('Evidence'
, 'Prior Probability for Theta'
, 'Posterior Probability for Theta'
, 'Most Plausible Theta'
), loc = 'upper left')
ax.set_xlabel("Theta")
ax.grid()
ax.set_title("Prior Distribution Updated with Some Evidence")
plt.show()

Now that we have a full distribution for the probability of various values of θ, we can take the mean of the distribution as our most plausible value for θ, which is about 0.27.

The data has caused us to believe that the true click-through rate is higher than we originally thought, but far lower than the 0.7 click-through rate observed so far from the facebook-yellow-dress campaign. Why is this the case? Note how wide our likelihood function is; it's telling us that there is a wide range of values of θ under which our data is likely. If the range of values under which the data were plausible were narrower, then our posterior would have shifted further. See what happens to the posterior if we observed a 0.7 click-through rate from 10, 100, 1,000, and 10,000 impressions:

import pymc3 as pm
import numpy as np

#create our data:
traces = {}
for ad_impressions in [10, 100, 1000, 10000]: #maintaining observed CTR of 0.7
clicks = np.array([ctr * ad_impressions])    #re-estimate the posterior for
impressions = np.array([ad_impressions])    #increasing numbers of impressions
with pm.Model() as model:
theta_prior = pm.Beta('prior', 11.5, 48.5)
observations = pm.Binomial('obs',n = impressions
, p = theta_prior
, observed = clicks)
start = pm.find_MAP()
step = pm.NUTS(state=start)
trace = pm.sample(5000
, step
, start=start
, progressbar=True)

f, ax = plt.subplots(1)
ax.plot(bins[:-1],prior_counts, alpha = .2)

counts = {}
for ad_impressions in [10, 100, 1000, 10000]:
posterior_counts, posterior_bins = np.histogram(trace['prior'], bins=[j/100. for j in xrange(100)])
posterior_counts = posterior_counts / float(len(trace))
ax.plot(bins[:-1], posterior_counts)
line0, line1, line2, line3, line4 = ax.lines
ax.legend((line0, line1, line2, line3, line4), ('Prior Distribution'
,'Posterior after 10 Impressions'
, 'Posterior after 100 Impressions'
, 'Posterior after 1000 Impressions'
,'Posterior after 10000 Impressions'))
ax.set_xlabel("Theta")
ax.axvline(ctr, linestyle = "--", alpha = .5)
ax.grid()
ax.set_ylabel("Probability of Theta")
ax.set_title("Posterior Shifts as Weight of Evidence Increases")
plt.show() 

As we obtain more and more data, we are more certain that the 0.7 success rate is the true success rate. Conditioning on more data as we update our prior, the likelihood function begins to play a larger role in our ultimate assessment because the weight of the evidence gets stronger. This would be particularly useful in practice if we wanted a continuous, fair assessment of how our campaigns are performing without having to worry about overfitting to a small sample.

There are a lot of concepts are beyond the scope of this tutorial, but are important for doing Bayesian analysis successfully, such as how to choose a prior, which sampling algorithm to choose, determining if the sampler is giving us good samplers, or checking for sampler convergence. Hopefully this tutorial inspires you to continue exploring the fascinating world of Bayesian inference.

Books:

Software:

Papers:

Other:

Want to keep learning? Download our new study from Forrester about the tools and practices keeping companies on the forefront of data science.

Photo by mattbuck. [CC BY-SA 3.0], via Wikimedia Commons

Aaron Kramer

Data analyst at DataScience. I'm into basketball, python, machine learning, algorithms and economics.