Probabilistic Programming and Bayesian Methods for Hackers using PyMC3 and Jupyter Notebooks¶

Getting Started¶

The easiest way to get started with an environment for Juypter Notebook, PyMC3, python (in version 3.6) and R is to install the anaconda package manager.

You will need to install several python packages as you can see from the import line above. You have two choices here, either you use anaconda or PyPi, the python package manager (pip install).

I personally use anaconda whenever possible and use PyPi if the package does not exist in anaconda. You should configure your anaconda environment to include the conda-forge anaconda channel.

I will not show it here, but it is possible to use R directly from within python via the rpy2 package. This is often quite handy if you do not have to leave the python environment.

In general you have 3 choices for doing "probabilistic programming" in python:

• PyMC3
• JAGS / PyJAGS
• Stan / PyStan

All packages implement Markov Chain Monte Carlo algorithms.

JAGS implements the so called Metropolis–Hastings algorithm. This can be used for all probability distributions, both, discrete and continuous, but will break down in high dimensional spaces.

Stan implements the so called Hamiltonian Monte Carlo (HMC) algorithm, where NUTS (The No U-Turn Sampler) is an optimized version of HMC and there is bascially never a reason to use plain HMC if you can use NUTS. HMC methods have the disadvantage that they only work for continuous probability distributions. But often you can find a way around this limitation. HMC/NUTS shines in high dimensional probability spaces and is typically "faster" than Metropolis-Hastings.

PyMC3 implements both algorithms above and allows for combining the different kinds of samplers.

All MCMC algorithms are sequential, e.g. it does not help you if you have many cores. But in order to verify if your MCMC results "converge" for correctness reasons you should ALWAYS run several chains in parallel and then use the $\hat R$ statistic (aka Gelman Rubin Diagnostic; should be close to 1.0) to verify the convergence. Michael Betancourt from the Stan team explains this here: Video: "Everything You Should Have Learned About Markov Chain Monte Carlo".

Bayes' Rule¶

Before we start we should have at least once explicitly stated Bayes' Rule. The key insight is this:

$$p(\Theta, D) = p(D, \Theta) = p(\Theta \;|\; D)\cdot p(D) = p(D \;|\; \Theta)\cdot p(\Theta)$$

Let's take this step by step.

Typically in the literature $\Theta$ stand for the set or group (this means that $\Theta$ typically does not stand for a single parameter, but for a group of parameters, e.g. $\{\gamma_1, \gamma_2, ..., \gamma_n\}$) of parameters of the model and $D$ stands for the observed data. What you have is the observed data. What you want to get to are the parameters that lead to these data observations.

The first two statements are just the joint probability distribution of parameters and data and because a probability distribution does not depend on the order of the parameters we can write $p(\Theta, D) = p(D, \Theta)$. That's basically the same joint distribution.

If you then use the definition of conditional probability distributions: $p(a, b) = p(a\;|\;b)\cdot p(b)$, which means that the joint distribution is the product of the "marginal distribution" $p(b)$ times the conditional distribution $p(a\;|\;b)$ conditional on $b$. That's just using the definition and nothing else.

Then you end up with Bayes' Rule:

$$\begin{align} p(\Theta \;|\; D)&= p(D \;|\; \Theta)\cdot p(\Theta) / p(D)\\ p(\Theta \;|\; D)&\propto p(D \;|\; \Theta)\cdot p(\Theta)\\ \end{align}$$

The second line means that $p(D)$ is typically not "needed". It is just a normalization constant to make sure that the probability distribution $p(D \;|\; \Theta)\cdot p(\Theta)$ is normalized to $1$. This also means that you can calculate $p(D)$ from $p(D \;|\; \Theta)\cdot p(\Theta)$ via:

$$p(D) = \int p(D \;|\; \Theta)\cdot p(\Theta) \hbox{d}\Theta$$ So all you really need is $p(D \;|\; \Theta)\cdot p(\Theta)$, which is also your model.

In words Bays' rule can be summarized as:

The probability distribution of the parameters of the model can be calculate by acting as if you know the parameters of the model and calculate the probability of the observed data item.

The probability of $\Theta$ given Data is proportional to the probability of Data given $\Theta$.

Warmup: Solving the Bayesian German Tank problem with PyMC3¶

The example is taken from here.

Here we will cover a classic problem in statistics, estimating the total number of tanks from a small sample. Suppose four tanks are captured with the serial numbers 10, 256, 202, and 97. Assuming that each tank is numbered in sequence as they are built, how many tanks are there in total?

The problem is named after the Allies' use of applied statistics in World War II to analyze German production capacity for Panzer tanks. The actual story of how that turned out is fascinating:

After the war, the allies captured German production records, showing that the true number of tanks produced in those three years was [...] almost exactly what the statisticians had calculated, and less than one fifth of what standard intelligence had thought likely.

Problem set-up¶

We start with a simplified story about how the data is generated — the generative model — and then attempt to infer specifics about that story (the parameters).

Here's our model:

1. Choose a number of tanks N to exist in the world.
2. Of those N, choose four that we happened to observe and call that data D.

The key is that some choices for N are much more likely than others. For example, if N was actually 5,000 then it would be a really odd coincidence that we happened to see all small numbers.

Solving it¶

Here's the gist: we want to know the posterior probability of the number of tanks N, given data D on which tank serial numbers have been observed so far:

$$P(N\;|\;D)\propto P(D\;∣\;N)\cdot P(N)$$

The right-hand side breaks down into two parts. The likelihood of observing all the serial numbers we saw is the product of the individual probabilities $D_i$ given the actual number of tanks in existence $N$.

$$\begin{align} P(D\;∣\;N)&=\Pi_i P(D_i\;∣\;N)\\ P(D_i\;∣\;N)&\sim \hbox{DiscreteUniform}(D, \hbox{min}=0, \hbox{max}=N)\\ \end{align}$$

The prior over N, which we know must be at least as high as the highest serial number we observed (call that $m$) but could be much higher:

$$P(N)\sim \hbox{DiscreteUniform}(N, \hbox{min}=m, \hbox{max}=\hbox{some big number})$$

PyMC3¶

So the mean here was 383, and the median was 323.

Linear Regression and Robust Linear Regression and the GLM¶

For the following examples I've combined content from different sources and tutorials on the web. Some of the sources I used were:

• PyMC3 documentation: linear regression
• PyMC3 documentation: robust linear regression
• Connect The Dots. Least Squares, Linear Regression, and Bayesian Regression

Theory of a basic Linear Regression Model¶

$$\vec y=\bf X\vec\beta+\vec\epsilon$$

With $\epsilon_i\sim \mathcal{N}(0,\sigma^2)$ = independent Gaussian noise on the measurement error.

For ordinary least squares (OLS) regression, we can quite conveniently solve this (find the $\beta$ values) using the Maximum Likelihood Estimate (MLE) which has a closed analytical form:

$$\hat\beta=({\bf X}^T{\bf X})^{-1}{\bf X}^T\vec y$$

Let's generate some toy data to play with:

So we generate data with an intercept of $a=5$ and an inclination of $b=2$.

Doing it with a pre-packaged scikit-learn solver¶

Doing it with PyMC3¶

Probabilistic Reformulation¶

Rather than using the formulation above:

$$\vec y=\bf X\vec\beta+\vec\epsilon$$ the probabilistic reformulation looks as follows: $$\vec y\sim\mathcal{N}(\bf X\vec\beta,\,\sigma\cdot {\bf I})$$

Robust Linear Regression¶

The only thing we've changed is in the glm_normal_likelihood = pm.Normal $\rightarrow$ glm_outliers_likelihood = pm.StudentT line. We've exchanged the normal distribution ($\mathcal{N}$) with the Student's T distribution. The normal distribution has "thin tails", while the Student's T distribution can have "fat tails" depending on the parameter $\nu$. The larger the parameter $\nu$ the closer the Student's T distribution is to the normal distribution. So a Student's T distribution with a small $\nu$ can "support" outliers much better than a normal distribution.

Please notice that the line that you perceive as a single blue line is actually a set of lines like in the graph above. Only because the $y$ scale is this time much larger you cannot visually se the difference that well. But if you could zoom in into the above graph you would see it.

Mixture Models¶

Finite Mixture Models¶

The following example is taken from the PyMC3 documentation.

Please note that on the right hand side of the above graph you can nicely see the effect of "burn-in". It takes roughly 2000 samples until the values "converge". This is why you typically throw away the first 2000 samples.

Infinite Mixture Models: Clustering the McDonalds Menu¶

This example is taken from here: Infinite Mixture Models with Nonparametric Bayes and the Dirichlet Process. This is also by far the best explanation I've found for the so called Stick-Breaking Process, and how it is equivalent to the Chinese Restaurant Process, the Polya Urn Model and the Dirichlet Process.

Have a look at Infinite Mixture Models with Nonparametric Bayes and the Dirichlet Process and scroll down to "Let’s dive into one of these clusterings" to see some of the clusters.

Have a look at scikit-learn's Gaussian mixture models documentation to learn more about this functionality.

Switchpoint analysis¶

Single switchpoint¶

The example is taken directly from the PyMC3 documentation: Case study 2: Coal mining disasters. You can look there for more details and explanations.

Multiple switchpoints: Bayesian Blocks / Data Binning / Data-Bucketing / Discretization¶

In a similar way like with our mixture model above where we treated the case with a fixed number of mixtures within PyMC3, but used pre-packaged algorithms for the infinite case, here I simply point you into the right direction.

We will look at a very specific case here: often you will need to discretize a continuous feature, also called binning or bucketing, e.g. you have some continuous data, but you want to make the representation more compact by treating the data as some countable number of cases.

There have been many rules of thumb used in the literature, e.g. Scott’s Rule and the Freedman-Diaconis Rule. Both are based on assumptions about the form of the underlying data distribution. A slightly more rigorous approach is Knuth’s rule, which is based on optimization of a Bayesian fitness function across fixed-width bins.

A further generalization of this approach is that of Bayesian Blocks, which optimizes a fitness function across an arbitrary configuration of bins, which need not be of equal width. This paper is quite accessible and the method can be understood as a switchpoint analysis with an infinite number of possible switchpoints.

The astroML package includes an implementation in python for all the above algorithms and an explanation of the algorithm is given in the blog post Dynamic Programming in Python: Bayesian Blocks.

Here is an example:

Supervised Learning: Data Binning / Data-Bucketing / Discretization¶

Just as a side remark, in case you have a supervised learning task where you want to classify data according to some features and you want to discretize some of the continuous features then a good approach for finding the bin boundaries is binning data using decision tree regression. The decision tree internally already builds up decision tree boundaries and those you can use as good bin boundaries points.

More theory is given in A Simple Guide to Entropy-Based Discretization.

Also have a look at decision trees in python with scikit-learn and pandas and its associated gist to learn more about how to use decision trees in scikit-learn.

Stochastic Volatility Model¶

The following example is taken from Probabilistic Programming in Python using PyMC by John Salvatier, Thomas Wiecki, Christopher Fonnesbeck. A similar example is also part of the PyMC3 documentation. And here is another similar discussion of this example.

One elegant example of how Probabilistic Programming can infer unobservable quantities of the stock market is the stochastic volatility model. Volatility is an important concept of quantitative finance as it relates to risk. Unfortunately, as Tony Cooper reminds us:

"Volatility is a strange thing - it exists but you can't measure it."

-- Tony Cooper on [Quantopian forums)(https://www.quantopian.com/posts/system-based-on-easy-volatility-investing-by-tony-cooper-at-double-digit-numerics).

If we can't measure something, the next best thing we can do is to try and model it. One way to do this is in a probabilistic framework is the concept of stochastic volatility: If we assume that returns are normally distributed, the volatility would be captured as the standard deviation of the normal distribution. Thus, the standard deviation gives rise to stochastic volatility. Intuitively we would assume that the standard deviation is high during times of market turmoil like the 2008 crash.

So the trivial thing to do would be to look at the rolling standard deviation of returns. But this is unsatisfying for multiple reasons:

• it often lags behind volatility,
• is a rather unstable measure,
• strongly dependent on the window size with no principled way of choosing it.

These properties are outlined in the plot below. The larger the window size, the more lag, the smaller the window size, the more unstable the estimate gets.

As we will see, the stochastic volatility model does a better job at all of them. But before we look at that we need to establish one more insight into the nature of volatility: it tends to cluster. This is nicely demonstrated by looking at the returns of the S&P 500 below. As you can see, during the 2008 financial crisis there is a lot volatility in the stock market (huge positive and negative daily returns) that gradually decreases over time.

So how do we model this clustering property? The Stochastic Volatility model assumes that that the standard-deviation of the returns follow a random-walk process. You can read the Wikipedia article I linked but essentially this process allows for slow, gradual changes over time.

What is interesting is that we can model the standard deviation itself to follow a random walk. Intuitively, we allow standard deviation to change over time but only ever so slightly at each time-point.

The plot further down below shows the latent volatility inferred from the model based on the market data. The lines represent the standard deviation of the Normal distribution we assume for the daily returns.

As you can see, the estimated volatility is much more robust. Moreover, all parameters are estimated from the data while we would have to find a reasonable window length for the rolling standard deviation ourselves. Finally, we do not just get a single estimate as with the rolling standard deviation but rather many solutions that are likely. This provides us with a measure of the uncertainty in our estimates and is represented by the width of the line below.

In summary, Bayesian statistics and Probabilistic Programming is an extremely powerful framework for Quantitative Finance as it provides:

• a principled framework for modeling latent, unobservable causes, as demonstrated by the stochastic volatility model;
• a measure of uncertainty in our estimates;
• powerful sampling methods that allow automatic estimation of highly complex models.

This example has 400+ parameters so using common sampling algorithms like Metropolis-Hastings would get bogged down, generating highly autocorrelated samples. Instead, we use NUTS, which is dramatically more efficient.

The Model¶

Asset prices have time-varying volatility (variance of day over day returns). In some periods, returns are highly variable, while in others they are very stable. Stochastic volatility models address this with a latent volatility variable, which changes over time. The following model is similar to the one described in the NUTS paper (Hoffman & Gelman, 2014, p. 21).

$$\begin{align} \sigma&\sim\hbox{exp}(50)\\ \nu&\sim\hbox{exp}(.1)\\ s_i&\sim\mathcal{N}(s_{i-1},\sigma)\\ log(\frac{y_i}{y_{i-1}})&\sim t(\nu, 0, e^{s_i})\\ \end{align}$$

Here, $y$ is the daily return series which is modeled with a Student-t distribution with an unknown degrees of freedom parameter, and a scale parameter determined by a latent process $s$. The individual $s_i$ are the individual daily log volatilities in the latent log volatility process.

Sometimes some of the probability distributions are parametrized with different parameters, e.g. in the below pm.GaussianRandomWalk may be parametrized in terms of $\tau$ where $\tau=1/\sigma^2$ (there is no documentation in the PyMC3 documenation about the pm.GaussianRandomWalk; you have to look this up in the sourcecode).

The same is true with different parametrizations of the Student T distribution, e.g. $\lambda=1/\sigma^2$.

Build the model¶

Our data consist of daily returns of the S&P 500 during the 2008 financial crisis.

Finite mixture model with time-varying probabilities¶

• Finite mixture model with time-varying probabilities

Bayesian Deep Learning¶

• Neural Networks in PyMC3 estimated with Variational Inference

Appendix: Further Reading¶

Usage of Bayesian Methods¶

By far the most instructive book that teaches you Bayesian methods is the book Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan 2nd Edition by John Kruschke. I highly recommend the book!

• Book Site with source.
• Book Blog

Some people have translated the examples in the book from R and JAGS to python and PyMC3:

• Doing Bayesian Data Analysis in Python
  • Osvaldo Martin: aloctavodia/Doing_bayesian_data_analysis: GitHub
    • jupyter nbviewer

I also can highly recommend to read the Stan manual from cover to cover. I postponed reading it for too long, because I thought it is a "boring" manual, but actually it explains quite a lot of modeling techniques!

Another quite good free online book is Probabilistic Programming & Bayesian Methods for Hackers by Cameron Davidson-Pilon (also available as a printed book)

While the book Machine Learning: A Probabilistic Perspective is not directly related to Markov Chain Monte Carlo methods it gives a good overview over many different probability models. The art in making best use of MCMC is in finding the right model and this book looks in detail into many different probability/likelihood models and especially important sets them into context and compares them. The book is not always easy to read and sometimes the notation is confusing. Nevertheless it helps a lot by pointing you into the right direction!

MCMC internals: how it all works¶

If you want to understand what is going on under the hood then I recommend to start with the following video from Michael Betancourt, a member of the Stan team: Video: "Everything You Should Have Learned About Markov Chain Monte Carlo".

And the paper A Conceptual Introduction to Hamiltonian Monte Carlo also by Michael Betancourt goes into more details. I suggest reading this paper together with Logic, Probability, and Bayesian Inference: A Conceptual Introduction to the Foundations of Applied Inference for Scientists and Engineers .

People to follow¶

The following list mentions some people who are very active in the MCMC area. It is worthwhile to look at some of the things they've done. Below I've given links to one of their web-sites, but often they are active in many different places in the web, e.g. simply use Google to search for them and their work:

• Michael Betancourt (some of his case studies are here)
• Austin Rochford
• Chris Fonnesbeck
• Thomas Wiecki (you can support him on patreon, so that he can do more quality work on PyMC3)
• Jim Savage (he is currently in the process of writing a book)
• Cam Davidson-Pilon

News from PyMC3¶

If you want ot learn more about PyMC3 then the best way is to follow some of the examples that others have done. Very noteworthy are the recent improvements in PyMC3 verion 3.1.

All around MCMC is an extremely active area of development. It is worthwhile to check regularly.