CS109B Data Science 2: Advanced Topics in Data Science

Lab 4 - Bayesian Analysis - EXTENDED¶

Harvard University
Spring 2020
Instructors: Mark Glickman, Pavlos Protopapas, and Chris Tanner
Lab Instructors: Chris Tanner and Eleni Angelaki Kaxiras
Content: Eleni Angelaki Kaxiras

---

Additions in this EXTENDED version¶

• Cleaner code in Bayesian Linear Regression. <---- Edited material
• Cleaner code in defining the Model <--- Edited Material
• Hierarchical Models. <----- New material
• Some info on Bayesian Logistic Regression with PyMC3 <-- New material

Learning Objectives¶

By the end of this lab, you should be able to:

• Apply Bayes Rule in calculating probabilities.
• Practice how to apply Bayesian analysis using PyMC3
• Avoid getting fired when talking to your Bayesian employer.

This lab corresponds to Lectures 6, 7, and 8, and maps to Homework 3.

Table of Contents¶

1. The Bayesian Way of Thinking or Is this a Fair Coin?
2. Intro to pyMC3.
3. Bayesian Linear Regression.

1. The Bayesian way of Thinking¶

Here is my state of knowledge about the situation. Here is some data, I am now going to revise my state of knowledge.

A. Bayes Rule¶

\begin{equation} \label{eq:bayes} P(A|\textbf{B}) = \frac{P(\textbf{B} |A) P(A) }{P(\textbf{B})} \end{equation}

$P(A|\textbf{B})$ is the posterior distribution, prob(hypothesis | data)

$P(\textbf{B} |A)$ is the likelihood function, how probable is my data B for different values of the parameters

$P(A)$ is the marginal probability to observe the data, called the prior, this captures our belief about the data before observing it.

$P(\textbf{B})$ is the marginal distribution (sometimes called marginal likelihood)

The problem we are about to solve gained fame as part of a game show "Let's Make A Deal" hosted by Monty Hall, hence its name. It was first raised by Steve Selvin in American Statistician in 1975.

The game is as follows: there are 3 doors behind one of which are the keys to a new car. There is a goat behind each of the other two doors. Let's assume your goal is to get the car and not a goat.

You are asked to pick one door, and let's say you pick Door1. The host knows where the keys are. Of the two remaining closed doors, he will always open the door that has a goat behind it. He'll say "I will do you a favor and open Door2". So he opens Door2 inside which there is, of course, a goat. He now asks you, do you want to open the initial Door you chose or change to Door3? Generally, in this game, when you are presented with this choice should you swap the doors? Initial Steps:

• Start by defining the events of this probabilities game. One definition is:

  • $A_i$: car is behind door $i$

  • $B_i$ host opens door $i$

$i\in[1,2,3]$

• In more math terms, the question is: is the probability of the price is behind Door 1 higher than the probability of the price is behind Door2, given that an event has occured?

B. Bayes Rule Revisited¶

We have data that we believe come from an underlying distribution of unknown parameters. If we find those parameters, we know everything about the process that generated this data and we can make inferences (create new data).

\begin{equation} \label{eq:bayes} P(\theta|\textbf{D}) = \frac{P(\textbf{D} |\theta) P(\theta) }{P(\textbf{D})} \end{equation}

But what is $\theta \;$?¶

$\theta$ is an unknown yet fixed set of parameters. In Bayesian inference we express our belief about what $\theta$ might be and instead of trying to guess $\theta$ exactly, we look for its probability distribution. What that means is that we are looking for the parameters of that distribution. For example, for a Poisson distribution our $\theta$ is only $\lambda$. In a normal distribution, our $\theta$ is often just $\mu$ and $\sigma$.

C. A review of Common Probability Distributions¶

Discrete Distributions¶

The random variable has a probability mass function (pmf) which measures the probability that our random variable will take a specific value $y$, denoted $P(Y=y)$.

• Bernoulli (binary outcome, success has probability $\theta$, $one$ trial): $ P(Y=k) = \theta^k(1-\theta)^{1-k} $

  ---

• Binomial (binary outcome, success has probability $\theta$, $k$ sucesses, $n$ trials): \begin{equation} P(Y=k) = {{n}\choose{k}} \cdot \theta^k(1-\theta)^{n-k} \end{equation}

Note: Binomial(1,$p$) = Bernouli($p$)

---

• Negative Binomial

  ---

• Poisson (counts independent events occurring at a rate $\lambda$) \begin{equation} P\left( Y=y|\lambda \right) = \frac{{e^{ - \lambda } \lambda ^y }}{{y!}} \end{equation} y = 0,1,2,...

  ---

• Discrete Uniform

  ---

• Categorical, or Multinulli (random variables can take any of K possible categories, each having its own probability; this is a generalization of the Bernoulli distribution for a discrete variable with more than two possible outcomes, such as the roll of a die)

  ---

• Dirichlet-multinomial (a generalization of the beta distribution for many variables)

Continuous Distributions¶

The random variable has a probability density function (pdf).

• Uniform (variable equally likely to be near each value in interval $(a,b)$) \begin{equation} P(X = x) = \frac{1}{b - a} \end{equation} anywhere within the interval $(a, b)$, and zero elsewhere.

  ---

• Normal (a.k.a. Gaussian) \begin{equation} X \sim \mathcal{N}(\mu,\,\sigma^{2}) \end{equation}

  A Normal distribution can be parameterized either in terms of precision $\tau$ or variance $\sigma^{2}$. The link between the two is given by \begin{equation} \tau = \frac{1}{\sigma^{2}} \end{equation}

  • Expected mean $\mu$
  • Variance $\frac{1}{\tau}$ or $\sigma^{2}$
  • Parameters: mu: float, sigma: float or tau: float

    ---

• Beta (where the variable ($\theta$) takes on values in the interval $[0,1]$, and is parametrized by two positive parameters, $\alpha$ and $\beta$ that control the shape of the distribution. Note that Beta is a good distribution to use for priors (beliefs) because its range is $[0,1]$ which is the natural range for a probability and because we can model a wide range of functions by changing the $\alpha$ and $\beta$ parameters.

\begin{equation} \label{eq:beta} P(\theta|a,b) = \frac{1}{B(\alpha, \beta)} {\theta}^{\alpha - 1} (1 - \theta)^{\beta - 1} \propto {\theta}^{\alpha - 1} (1 - \theta)^{\beta - 1} \end{equation}

where the normalisation constant, $B$, is a beta function of $\alpha$ and $\beta$,

\begin{equation} B(\alpha, \beta) = \int_{t=0}^1 t^{\alpha - 1} (1 - t)^{\beta - 1} dt. \end{equation}

---

• Exponential

  ---

• Gamma

Code Resources:¶

• Statistical Distributions in numpy/scipy: scipy.stats
• Statistical Distributions in pyMC3: distributions in PyMC3 (we will see those below).

Poisson¶

Change the value of $\lambda$ in the Poisson PMF and see how the plot changes. Remember that the y-axis in a discrete probability distribution shows the probability of the random variable having a specific value in the x-axis.

\begin{equation} P\left( X=y \right|\lambda) = \frac{{e^{ - \lambda } \lambda ^y }}{{y!}} \end{equation}

for $y \ge0$.

Routine is stats.poisson.pmf(x, lambda). $\lambda$ is our $\theta$ in this case. $\lambda$ is also the mean in this distribution.

Binomial¶

Uniform¶

Change the value of $\mu$ in the Uniform PDF and see how the plot changes.

Remember that the y-axis in a continuous probability distribution does not shows the actual probability of the random variable having a specific value in the x-axis because that probability is zero!. Instead, to see the probability that the variable is within a small margin we look at the integral below the curve of the PDF.

The uniform is often used as a noninformative prior.

Uniform - numpy.random.uniform(a=0.0, b=1.0, size)

$\alpha$ and $\beta$ are our parameters. size is how many tries to perform. Our $\theta$ is basically the combination of the parameters a,b. We can also call it \begin{equation} \mu = (a+b)/2 \end{equation}

Beta¶

If we apply the formula we can see why we get these plots for the various values of $a$ and $b$.
Notice that for $a=b=1.$ we get a constant. From then on, as the values increase, we get a curve that looks more and more like a Gaussian.

Gaussian¶

D. Is this a Fair Coin?¶

Let's say you visit the casino in Monte Carlo. You want to test your theory that casinos are dubious places where coins have been manipulated to have a larger probability for tails. So you will try to estimate how fair a coin is based on a certain amount of flips.
You have no prior opinion on the coin's fairness (i.e. what $p$ might be), and begin flipping the coin. You get either Heads ($H$) or Tails ($T$) as our observed data and want to see if your posterior probabilities change as you obtain more data, that is, more coin flips. A nice way to visualize this is to plot the posterior probabilities as we observe more flips (data).

We will be using Bayes rule. $\textbf{D}$ is our data.

\begin{equation} \label{eq:bayes} P(\theta|\textbf{D}) = \frac{P(\textbf{D} |\theta) P(\theta) }{P(\textbf{D})} \end{equation}

In the case of a coin toss when we observe $k$ heads in $n$ tosses:

\begin{equation} \label{eq:bayes} P(\theta|\textbf{k}) = Beta(\alpha + \textbf{k}, \beta + n - \textbf{k}) \end{equation}

we can say that $\alpha$ and $\beta$ play the roles of a "prior number of heads" and "prior number of tails".

Top

2. Introduction to pyMC3¶

PyMC3 is a Python library for programming Bayesian analysis, and more specifically, data creation, model definition, model fitting, and posterior analysis. It uses the concept of a model which contains assigned parametric statistical distributions to unknown quantities in the model. Within models we define random variables and their distributions. A distribution requires at least a name argument, and other parameters that define it. You may also use the logp() method in the model to build the model log-likelihood function. We define and fit the model.

PyMC3 includes a comprehensive set of pre-defined statistical distributions that can be used as model building blocks. Although they are not meant to be used outside of a model, you can invoke them by using the prefix pm, as in pm.Normal.

Markov Chain Monte Carlo (MCMC) Simulations¶

PyMC3 uses the No-U-Turn Sampler (NUTS) and the Random Walk Metropolis, two Markov chain Monte Carlo (MCMC) algorithms for sampling in posterior space. Monte Carlo gets into the name because when we sample in posterior space, we choose our next move via a pseudo-random process. NUTS is a sophisticated algorithm that can handle a large number of unknown (albeit continuous) variables.

Top

3. Bayesian Linear Regression¶

Defining the Problem¶

Our problem is the following: we want to perform multiple linear regression to predict an outcome variable $Y$ which depends on variables $\bf{x}_1$ and $\bf{x}_2$.

We will model $Y$ as normally distributed observations with an expected value $mu$ that is a linear function of the two predictor variables, $\bf{x}_1$ and $\bf{x}_2$.

\begin{equation} Y \sim \mathcal{N}(\mu,\,\sigma^{2}) \end{equation}

\begin{equation} \mu = \beta_0 + \beta_1 \bf{x}_1 + \beta_2 x_2 \end{equation}

where $\sigma^2$ represents the measurement error (in this example, we will use $\sigma^2 = 10$)

We also choose the parameters to have normal distributions with those parameters set by us.

\begin{eqnarray} \beta_i \sim \mathcal{N}(0,\,10) \\ \sigma^2 \sim |\mathcal{N}(0,\,10)| \end{eqnarray}

We will artificially create the data to predict on. We will then see if our model predicts them correctly.

Artificially creating some data to test our model.¶

Now let's see if our model will correctly predict the values for our unknown parameters, namely $b_0$, $b_1$, $b_2$ and $\sigma$.

Defining a Model in PyMC3¶

Building the model¶

Step1: We choose the probability model for our experiment.

Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

Step2: Choose a prior distributions of the two rates, what we believe the rates were before we observed the data, and the switchpoint.

 beta0 = Normal('beta0', mu=0, sd=10)
 # Note: betas is a vector of two variables, b1 and b2, (denoted by shape=2)
 # so, in array notation, our beta1 = betas[0], and beta2=betas[1]
 betas = Normal('betas', mu=0, sd=10, shape=2) 
 sigma = HalfNormal('sigma', sd=1)

Note: Watch for missing values. Missing values are handled transparently by passing a MaskedArray or a pandas.DataFrame. Behind the scenes, another random variable, disasters.missing_values is created to model the missing values. If you pass a np.array with missing values you will get an error.

Note: If our problem was a classification for which we would use Logistic regression see below

Fitting the Model with MAP (FYI, we will not directly use this method)¶

In Bayesian analysis we have our prior(s), we define our likelihood, and, having specified our model, we try to calculate posterior estimates for the unknown variables in the model. We could try to calculate the posterior estimates analytically, but for most the models, this is not feasible. What we do then is compute summaries based on samples drawn from the posterior distribution using Markov Chain Monte Carlo (MCMC) sampling methods.

\begin{equation} \label{eq:beta} P(\theta|\textbf{D}) \rightarrow \{\theta_1,....\theta_n\} \end{equation}

Then we can find any estimate we want by using these samples, for example:

\begin{equation} \mathbb{E}[f(\theta] = \int d\theta{p(\theta) f(\theta)} \end{equation}

So we calculate the maximum a posteriori (MAP) point using optimization methods.

\begin{equation} f(\hat{\theta}), \hat{\theta} = argmax ({p(\theta))} \end{equation}

The maximum a posteriori (MAP) estimate for a model, is the mode of the posterior distribution and is generally found using numerical optimization methods. This is often fast and easy to do, but only gives a point estimate for the parameters and can be biased if the mode isn’t representative of the distribution. PyMC3 provides this functionality with the find_MAP function.

MAP estimate is not always reasonable, especially if the mode is at an extreme or we have a multimodal distribution, or we have high dimensional posteriors. This will often occur in hierarchical models with the variance parameter for the random effect. If the individual group means are all the same, the posterior will have near infinite density if the variance parameter for the group means is almost zero. Most techniques for finding the MAP estimate only find a local optimium (which is often good enough), and can therefore fail badly for multimodal posteriors, as mentioned above.

To solve these issues we turn to sampling as our method for finding the posterior.

You do not have to worry about MAP in our problems. Our pyMC3 models use the MAP method to initialize the variables under the hood and we do not have to explicitly set this.

Fitting the Model with Sampling - Doing Inference¶

See below for PyMC3's sampling method. As you can see it has quite a few parameters. Most of them are set to default values by the package. For some, it's useful to set your own values.

pymc3.sampling.sample(draws=500, step=None, init='auto', n_init=200000, start=None, trace=None, chain_idx=0, chains=None, cores=None, tune=500, progressbar=True, model=None, random_seed=None, discard_tuned_samples=True, compute_convergence_checks=True, **kwargs)

Parameters to set:

• draws (int): number of samples to draw, defaults to 500.
• tune (int): number of iterations to tune, defaults to 500.
• target_accept (float in $[0, 1]$). The step size is tuned such that we approximate this acceptance rate. Higher values like 0.9 or 0.95 often work better for problematic posteriors.
• (optional) cores (int) number of chains to run in parallel, defaults to the number of CPUs in the system, but at most 4.

pm.sample returns a pymc3.backends.base.MultiTrace object that contains the samples. We usually name it trace. All the information about the posterior is in trace, which also provides statistics about the sampler.

Model Plotting¶

PyMC3 provides a variety of visualizations via plots: https://docs.pymc.io/api/plots.html.

One of them is the traceplot, another is the compareplot.

$\hat{R}$ is a metric for comparing how well a chain has converged to the equilibrium distribution by comparing its behavior to other randomly initialized Markov chains. Multiple chains initialized from different initial conditions should give similar results. If all chains converge to the same equilibrium, $\hat{R}$ will be 1. If the chains have not converged to a common distribution, $\hat{R}$ will be > 1.01. $\hat{R}$ is a necessary but not sufficient condition.

For details on the $\hat{R}$ see Gelman and Rubin (1992).

This linear regression example is from the original paper on PyMC3: Salvatier J, Wiecki TV, Fonnesbeck C. 2016. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2:e55 https://doi.org/10.7717/peerj-cs.55

Top

4. Hierarchical Models¶

Gelman et al. famous radon dataset is a classic for hierarchical modeling. In this dataset the amount of the radioactive gas radon has been measured among different households in all county's of several states. Radon gas is known to be the highest cause of lung cancer in non-smokers. It is believed to be more strongly present in households containing a basement and to differ in amount present among types of soil. Here we'll investigate this differences and try to make predictions of radonlevels in different county's based on the county itself and the presence of a basement.

Each row in the dataframe represents the radon measurements for one house in a specific county including whether the house has a basement (floor = 0) or not (floor = 1). We are interested in whether having a basement increases the radon measured in the house.

To keep things simple let's keep only the following three variables: county, log_radon, and floor

Let's check how many different counties we have. We also notice that they have a different number of houses. Some have a large number of houses measured, some only 1.

Pooling: Same Linear Regression for all¶

We can just pool all the data and estimate one big regression to asses the influence of having a basement on radon levels across all counties. Our model would be:

\begin{equation} y_{i} = \alpha + \beta*floor_{i} \end{equation}

Where $i$ represents the measurement (house), and floor contains a 0 or 1 if the house has a basement or not. By ignoring the county feature, we do not differenciate on counties.

Remember, with the pooled model we have only one intercept, $\alpha$, and only one slope, $\beta$ for all the counties. Let's plot the regression lines.

Unpooling: Separate Linear Regression for each county¶

We believe that different counties have different relationships of radon and basements. Our model would be:

\begin{equation} radon_{i,c} = \alpha_c + \beta_c*floor_{i,c} \end{equation}

Where $i$ represents the measurement, $c$ the county, and floor contains a 0 or 1 if the house has a basement or not.

Notice we have separate coefficients for each county in $a_c$ and $b_c$. They are totally different, they do not even come from the same distribution.

We will do this for only one county, as an example. We pick HENNEPIN county.

Print the regression line for our chosen county alone.

Partial pooling: Hierarchical Regression (Varying-Coefficients Model)¶

Merely by the fact that all counties are counties, they share similarities, so there is a middle ground to both of these extremes. Specifically, we may assume that while $\alpha_c$ and $\beta_c$are different for each county as in the unpooled case, the coefficients are all drawn from the same distribution:

\begin{equation} radon_{i,c} = \alpha_c + \beta_c*floor_{i,c} \end{equation}

\begin{equation} a_c \sim \mathcal{N}(\mu_a,\,\sigma_a^{2}) \end{equation}

\begin{equation} b_c \sim \mathcal{N}(\mu_b,\,\sigma_b^{2}) \end{equation}

where the common parameters are: \begin{eqnarray} \mu_a \sim \mathcal{N}(0,\,10) \\ \sigma_a^2 \sim |\mathcal{N}(0,\,10)| \\ \mu_b \sim \mathcal{N}(0,\,10) \\ \sigma_b^2 \sim |\mathcal{N}(0,\,10)| \end{eqnarray}

Add math for mu and sigma

The different counties are effectively sharing information through the commin priors. We are thus observing what is known as shrinkage; modeling the groups not as independent from each other, neither as a single group but rather as related.

Defining the Model for the Hierarchical Model¶

Inference¶

This tutorial has borrowed from PyMC3 docs: https://docs.pymc.io/notebooks/GLM-hierarchical.html

What about Logistic Regression?¶

If the problem above was a classification that required a Logistic Regression, we would use the logistic function ( where $\beta_0$ is the intercept, and $\beta_i$ (i=1, 2, 3) determines the shape of the logistic function).

\begin{equation} Pr(Y=1|X_1,X_2,X3) = {\frac{1}{1 + exp^{-(\beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3)}}} \end{equation}

Since both $\beta_0$ and the $\beta_i$s can be any possitive or negative number, we can model them as gaussian random variables.

\begin{eqnarray} \beta_0 \sim \mathcal{N}(\mu,\,\sigma^2) \\ \beta_i \sim \mathcal{N}(\mu_i,\,\sigma_i^2) \end{eqnarray}

In PyMC3 we can model those as:

pm.Normal('beta_0', mu=0, sigma=100)

(where $\mu$ and $\sigma^2$ can have some initial values that we assign them, e.g. 0 and 100)

The dererministic variable would be:

p-logit = beta0 + beta_1 * X_1 + beta_2 * X_2 + beta_3 * X_3

To connect this variable (p_logit) with our observed data, we would use a Bernoulli as our likelihood.

our_likelihood = pm.Bernoulli('our_likelihood', logit_p=p-logit, observed=our_data)

Notice that the main difference with Linear Regression is the use of a Bernoulli distribution instead of a Gaussian distribution, and the use of the logistic function instead of the identity function.

Note: we could explicitly create this variable using pm.Deterministic.

    mu = pm.Deterministic('mu', beta0 + beta_1*X_1 + beta_2*X_2 + beta_3*X_3)

References:

• Salvatier J, Wiecki TV, Fonnesbeck C. 2016. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2:e55 (https://doi.org/10.7717/peerj-cs.55)
• Distributions in PyMC3
• More Details on Distributions

Information about PyMC3 functions including descriptions of distributions, sampling methods, and other functions, is available via the help command.

Cool Reading¶

• How Bayesian Analysis and Lawrence D. Stone found the Wreckage of Air France Flight AF 447.
• Search for the gold on the sunken SS Central America.