import numpy as np
import pandas as pd
from scipy import stats
from IPython.display import YouTubeVideo
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
Bayesian Inference#
In this section, we’ll focus on computing and using posterior distributions in more sophisticated Bayesian models. We’ll start by discussing why posterior distributions are useful in Bayesian inference, and then explain why they’re hard. Then, in the next section, we’ll learn about approximating distributions using sampling.
Why we need posterior distributions#
In general, we need the posterior distribution so that we can make statements and decisions about our unknown quantity of interest,
In the product review example:
Our parameter of interest
represents the probability of a positive review.If we chose a Beta prior, i.e.,
, then the posterior distribution also belonged to the Beta family: .This made it easy to determine things like the MAP estimate or LMSE estimate, simply by using known properties of the Beta distribution.
But what if our posterior distribution didn’t have such a convenient form? In that case, we would have to compute the posterior (and any estimates from it) ourselves:
\begin{align} p(\theta|x) &= \frac{p(x|\theta)p(\theta)}{p(x)} \ &= \frac{p(x|\theta)p(\theta)}{\int p(x|\theta)p(\theta),d\theta} \ \end{align}
In general, the integral in the denominator could be impossible to compute. We call the denominator the normalizing constant: it’s a constant because it doesn’t depend on
In the next section, we’ll see a few examples that illustrate why computing the normalizing constant is hard, but first, let’s look at three examples of why we might need to know it in the first place.
Computing probabilities#
Suppose we want to know the probability that
\begin{align} P(\theta > 0.7 | x) &= \int_{0.7}^1 p(\theta|x) , dx \ &= \int_{0.7}^1 \frac{p(x|\theta)p(\theta)}{p(x)} , dx \ &= \frac{1}{p(x)} \int_{0.7}^1 p(x|\theta)p(\theta) , dx \end{align}
In the last step, we used the fact that p(x) doesn’t depend on
If we don’t know
MAP Estimation#
Suppose we want to compute the MAP estimate:
In the last step, we used the fact that p(x) doesn’t depend on
If
To summarize: for low-dimensional continuous variables, or discrete random variables with a low number of possible values, we can compute the MAP estimate without needing to know the exact posterior. For higher-dimensional random variables and/or discrete random variables with many possible values, this won’t work.
LMSE Estimation#
Suppose we want to compute the LMSE estimate. Recall the definition of conditional expectation (see Data 140 textbook, Chapter 9 and Chapter 15):
In order to compute the LMSE estimate, we need to compute the denominator,
The same is true for computing the expected value of any other function of
To summarize: any computations involving the posteriors (probabilities, expectations, etc.) require us to have the full normalized distribution: the numerator in Bayes’ rule isn’t enough.
Why computing posterior distributions is hard#
In simple models like our product review model or our model for heights, it was easy to compute the exact posterior for the unknown variable that we were interested in. This happened because we chose a conjugate prior. In most other cases, computing the exact posterior is hard! Here are two examples:
One-dimensional non-conjugate prior#
Let’s return to the product review example, but this time, instead of a Beta prior, we choose
This distribution looks much more complicated: we can’t reduce it to a known distribution at all. So, in order to properly compute
This integral is difficult to solve in closed form. In this specific example, since this is a one-dimensional problem, we could take advantage of numerical integration. In other words, for a particular sequence of values
Multi-dimensional example#
Consider the exoplanet model from last time:
We can write the likelihood and prior. To simplify, we’ll write
We can try computing the posterior over the hidden variables
This distribution is more complicated than anything we’ve seen up until now. It’s the joint distribution over
Computing the normalization constant
For our dataset of over 500 planets, the sums alone would require a completely intractable amount of computation:
2**517
429049853758163107186368799942587076079339706258956588087153966199096448962353503257659977541340909686081019461967553627320124249982290238285876768194691072
Worse still, we can’t even compute the MAP estimate for the labels
Even in this fairly simple model, with two groups, we’ve found that exact inference is completely hopeless: there’s no way we can compute the exact posterior for all our unknowns. In the rest of this section, we’ll talk about ways to get around this problem using approximations to the posterior distribution.
Specifically, our approximation methods will take advantage of what we’ve learned. We know that the hardest part of computing posterior distributions is computing the normalization constant