Title :¶

Bayes PyMC3

Description :¶

Run this notebook in your own environment. Do not forget to download the images/ and data/ folders.

CS109B Data Science 2: Advanced Topics in Data Science

Lecture 10 - Bayesian Analysis and Introduction to pyMC3¶

Harvard University
Spring 2021
Instructors: Pavlos Protopapas, Mark Glickman, and Chris Tanner
Additional Instructor and content: Eleni Angelaki Kaxiras

---

Learning Objectives¶

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

• identify and describe some of the most important probability distributions.
• apply Bayes Rule in calculating probabilities (and would have had fun in the process).
• create probabilistic models in the PyMC3 library.

Table of Contents¶

1. Some common families of distributions (review)
2. The Bayesian Way of Thinking
3. Bayesian Regression with pyMC3 and Defining a model in pyMC3
4. Bayesian Logistic Regression with pyMC3
   
   Appendix

1. Some common families of distributions¶

Statistical distributions are characterized by one of more parameters, such as $\mu$ and $\sigma^2$ for a Gaussian distribution.

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

---

Discrete Distributions¶

The probability mass function (pmf) of a discrete random variable $Y$ is equal to the probability that our random variable will take a specific value $y$: $f_y=P(Y=y)$ for all $y$. The variable, most of the times, assumes integer values. Plots for pmfs are better understood as stem plots since we have discrete values for the probabilities and not a curve. Probabilities should add up to 1 for proper distributions.

• Bernoulli for a binary outcome, success has probability $\theta$, and we only have $one$ trial:
  

\begin{equation} P(Y=k) = \theta^k(1-\theta)^{1-k} \end{equation}

• Binomial for a binary outcome, success has probability $\theta$, $k$ successes in $n$ trials:

Our random variable is $Y$= number of successes in $n$ trials.

\begin{equation} P(Y=k|n,\theta) = {{ n}\choose{k}} \cdot \theta^k(1-\theta)^{n-k} \quad k=0,1,2,...,n \end{equation}

As a reminder ${{ n}\choose{k}}$ is "from $n$ choose $k$": \begin{equation} {{ n}\choose{k}} = \frac{n!}{k!(n-k)!} \end{equation}

$EY=n\theta$, $VarY = np(1-p)$

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

• Negative Binomial for a binary outcome, success has probability $\theta$, we have $r$ successes in $x$ trials:

Our random variable is $X$= number of trials to get to $r$ successes.

\begin{equation} P(X=x|r,\theta) = {{ x-1}\choose{r-1}} \cdot \theta^r(1-\theta)^{x-r} \quad x=r,r+1,... \end{equation}

• Poisson (counts independent events and has a single parameter $\lambda$)
  

\begin{equation} P\left( Y=y|\lambda \right) = \frac{{ e^{ - \lambda } \lambda ^y }}{{ y!}} \quad y = 0,1,2,... \end{equation}

• Discrete Uniform for a random variable $X\in(1,N)$:
  

\begin{equation} P(X=x|N) = \frac{1}{N}, \quad x=1,2,...,N \end{equation}

• 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)

\begin{equation} f(x_1,x_2,...,x_n) = \frac{m}{x_1!\cdot x_2!\dotsb x_n!} \cdot p_1^{x_1}\cdot p_2^{x_2}\dotsb p_n^{x_n} \end{equation}

---

Continuous Distributions¶

The random variable has a probability density function (pdf) whose area under the curve equals to 1.

• Uniform (variable equally likely to be near each value in interval $(a,b)$) \begin{equation} f(x|a,b) = \frac{1}{b - a} \quad x\in [a,b] \quad \text{and 0 elsewhere}. \end{equation}

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. The mean is $\mu = \frac{(a+b)}{2}$

• 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 value (mean) $\mu$
  • Variance $\frac{1}{\tau}$ or $\sigma^{2}$
  • Parameters: mu: float, sigma: float or tau: float
  • Range of values (-$\infty$, $\infty$)

• 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. Its density (pdf) is:

\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_{x=0}^1 x^{\alpha - 1} (1 - x)^{\beta - 1} dx. \end{equation}

• 'Nice', unimodal distribution
• Range of values $[0, 1]$

$EX = \frac{a}{a+b}$, $VarX = \frac{ab}{(a+b)^2(a+b+1)}$

Code Resources:¶

• Statistical Distributions in numpy/scipy: scipy.stats

Top

2. 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.

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, probability(parameters| data).

$P(\textbf{B} |A)$ is the likelihood function, how probable is my data 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). This serves a normalization purpose.

Let's Make a Deal¶

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 Door1 you chose or change to Door3? Generally, in this game, when you are presented with this choice should you swap the doors?

Hint:

• 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?

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} 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$.

Top

3. Bayesian Regression with 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. They are not meant to be used outside of a model, and you can invoke them by using the prefix pm, as in pm.Normal.

For more see: PyMC3 Quickstart

Distributions in PyMC3:¶

• Statistical Distributions in pyMC3.

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

Defining a Model in PyMC3¶

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.

Let's create some synthetic data¶

Building the model¶

Step1: Formulate the probability model for our data: $Y \sim \mathcal{N}(\mu,\,\sigma^{2})$

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

Step2: Choose a prior distribution for our unknown parameters.

 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)

Step3: Construct the likelihood function.

Step4: Determine the posterior distribution, this is our main goal.

Step5: Summarize important features of the posterior and/or plot the parameters.

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

Now all we need to do is sample our model.

... to be continued¶

Appendix A: Bayesian Logistic Regression with pyMC3¶

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:

logitp = beta0 + beta_1 * X_1 + beta_2 * X_2 + beta_3 * X_3

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

our_likelihood = pm.Bernoulli('our_likelihood', logit_p=logitp, 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.

Appendix B: Is this a fair coin?¶

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} P(\theta|\textbf{D}) = \frac{P(\textbf{D} |\theta) P(\theta) }{P(\textbf{D})} \end{equation}

We start with a non-informative prior, a Beta distribution with (a=b=1.)

\begin{equation} P(\theta|\textbf{k=0}) = Beta(1., 1.) \end{equation}

Then, as we get new data (say, we observe $k$ heads in $n$ tosses), we update our Beta with new a,b as follows:

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

(the proof is beyond our scope, if interested, see this Wikipedia article)

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

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

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

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.