Simulation and Bayes Rule

Here, I'll briefly describe some of the step methods for producing a Monte Carlo Markov Chain. We have discussed in detail the random walk Metropolis Hastings sampler. We will briefly discuss in a bit more detail the following samplers:

1. Gibbs Samplers
2. Metropolis Hastings Samplers
   a. Random Walk Metropolis Hastings
   b. Adaptive Metropolis Hastings (we won't cover)
   c. Independent Metropolis Hastings (we won't cover)
   d. Slice Samplers and Other variants (we won't cover)

The nitty-gritty details of the following two samplers are beyond the scope of the course, but I will try to give you a feel for how they improve on Metropolis Hastings because they are known to be a major improvement over any of the samplers listed above even for high dimensional models (lots of coefficients to estimate):

• Hamiltonian Monte Carlo Sampler (No U-turn Sampler)
• Affine Invariant Sampler (EMCEE)

This is a very fast moving topic- so fast, in fact- that the current (non-beta) version of pymc doesn't have either NUTS or the Affine Invariant Sampler.

Gibbs Sampler¶

The Gibbs Sampler is a special case of the Random-Walk Metropolis Hastings algorithm and one worth knowing about in a bit of detail because many tutorials and discussions about MH (especially older ones) are entertwined with discussions on Gibbs sampling and it can be confusing for the uninitiated. Despite this being a special case of MH, it is worth seeing how it works in a bit of detail. The idea is to successively sample each $\theta$ from its conditional distribution:

\begin{align} \theta^{t+1}_1 \sim & f_1(\theta^{t+1}_1 \mid \theta^{t}_2,\theta^{t}_3,\ldots, \theta^{t}_k,\mathbf{y},\mathbf{x}) \\ \theta^{t+1}_2 \sim & f_2(\theta^{t+1}_2 \mid \theta^{t+1}_1,\theta^{t}_3,\ldots, \theta^{t}_k,\mathbf{y},\mathbf{x}) \\ \theta^{t+1}_3 \sim & f_3(\theta^{t+1}_3 \mid \theta^{t+1}_1,\theta^{t+1}_2,\ldots, \theta^{t}_k,\mathbf{y},\mathbf{x})\\ &\vdots \\ \theta^{t+1}_k \sim & f_k(\theta^{t+1}_k \mid \theta^{t+1}_1,\theta^{t+1}_2,\theta^{t+1}_3\ldots, \theta^{t+1}_{k-1},\mathbf{y},\mathbf{x}) \end{align}

where each parameter $i$ has a conditional distribution function $f_i$.

First consider this simple example, where our parameters to be estimate are $\theta = [\mu,\sigma]$ and where we have flat priors:

When we employ a Gibbs Sampler, all samples are accepted (there are versions of the Gibbs Sampler that use a rejection based approach- not covered here). This is accomplished by choosing problems where the conditional distributions:

\begin{align} f_\mu(\mu | \sigma, \mathbf{Y}) \\ f_\sigma(\sigma | \mu, \mathbf{Y}) \end{align}

are known and have closed form solutions. For example, if we have

\begin{equation} f_\mu(\mu | \sigma, \mathbf{Y}) = \phi \left ( \frac{\mathbf{Y} - \mu}{\sigma} \right ) \end{equation}

that is, Y is distributed normal with mean $\mu$ and standard deviation $\sigma$, then it can be shown that \begin{equation} f_\sigma \left(\sigma^2 | \mu, \mathbf{Y}\right ) = InverseGamma(k,s) \end{equation}

where $k$ and $s$, the shape and scale parameters for the Inverse Gamma Distribution are \begin{align} k &= \frac{n}{2} \\ s &= \left (\frac{(Y-\mu)^2}{2} \right )^{-1} \end{align}

The Gibbs Sampler then proceeds as follows:

1. Start the sampler at an arbitrary starting value for $\theta$. In my code below, I start the Gibbs sampler at $\theta_0 = [\mu_0=30,\sigma_0=3]$.
2. Draw the next (always accepted) value of $\sigma_t$ by:
   1. Draw random variate from a $\Gamma(k,s) = \Gamma \left ( \frac{n}{2},\left (\frac{(Y-\mu_{t-1})^2}{2} \right )^{-1} \right )$ and denote as $g$
   2. Define $var_t = \frac{1}{g}$
3. Draw the next (always accepted) value of $\mu_t$ by:
   1. Draw a normal random variate having mean $\bar{Y}$ and standard deviation $\sqrt{\frac{var_t}{n}}$ and denote as $\mu_t$.
4. Store the values $[\mu_t,\sqrt{var_t}]$ as the next sampled values for $\mu$ and $\sigma$.

Important thing to note about this process:

1. Gibbs sampling depends on our ability to analytically solve for conditional distributions.
2. This involves derivations in mathematical statistics
3. Once we have these closed form solutions, we can accept every sampled draw

The Gibbs Sampler was the predominant sampling method early on in applied Bayesian statistics because it was the empirical analogue of conjugate priors (the focus of Bayesian Statistics before the computer age) and does have real advantages over MH Random Walk for problems having solvable conditional distributions, since it accepts every sample and can be more efficient (less computing time to reach convergence and inference).

Despite this, both the Gibbs Sampler and the MH Sampler are less than ideal for

1. High dimensional problems (many parameters)
2. High correlation amongst parameters

Bayesian OLS Regression applied to Tobias and Koop¶

We now turn our attention back to MH Random Walk sampling using a real dataset.

This is the model statement describing priors and the likelihood. Here, we have 4 beta variables and sigma defines stochastic variables (we want a chain of sampled values for all of these variables) and we provide a prior distribution and hyper-parameters for each variable. The likelihood function is chosen to be Normal, having mean denoted as mu (which is defined as $\mathbf{x}\beta$, having standard deviation $\sigma$ (denoted as sigma). Our "dependent variable" is given by observed=tobias_koop.ln_wage, where tobias_koop.ln_wage is the log of wage from the Tobias and Koop dataset loaded above. In more formal terms, the code below sets up a basic_model having the following form:

\begin{align} Prob(\mu|\sigma,\mathbf{data}) \propto& Prob(\mathbf{data}|\mu,\sigma) \times Prob(\sigma | \mu^0_\sigma, \sigma^0_\sigma) \times Prob(\text{prior on }\beta) \\ Prob(\mu|\sigma,\mathbf{data}) \propto& \phi \left( \frac{\mathbf{ln\_wage} - \mathbf{x}\beta}{\sigma} \right ) \times HalfNormal \left(\frac{\sigma - \mu^0_\sigma}{\sigma^0_\mu} \right ) \times C \\ Prob(\mu|\sigma,\mathbf{data}) \propto& \phi \left( \frac{\mathbf{ln\_wage} - \mathbf{x}\beta}{\sigma} \right ) \times HalfNormal \left(\frac{\mu - 0}{10} \right ) \end{align}

where $\phi(.)$ is the standard normal pdf. The constant C (the joint probability of the prior on our $\beta$'s) is invariant with changes in our chain, and can therefore be dropped.

Metropolis-Hastings Random Walk¶

I like visualizing my distributions using the corner plot function. To do that, convert your chains into a pandas dataframe:

What is the problem? We have what appears to be a converged chain that gives us back what we expect if we compare to the OLS results. Diving into the details of the sampler, one can see that MH Random Walk, even for this relatively simple model has some pairs of parameters that have high covariance. Focusing on the Intercept and Educ, the following video shows the evolution of our chain for each of these two variables, and how it explores the probability surface:

Hamiltonian Monte Carlo and the No U-turn Sampler¶

Hamiltonian Monte Carlo was developed initially in the field of particle physics for simulating the movement of molecules in three dimensions. It is trivially extended to $k$ dimensions (our parameter space). This step method adds an auxiliary parameter for each model parameter and therefore doubles the parameter space in a way that actually makes the exploration of the posterior more efficient. These auxiliary parameters (nuisance parameters) are discarded after estimation. Our original model parameters (helps us to define potential energy which is the posterior probability) and these auxiliary parameters (helps us define kinetic energy), which aid in a more efficient exploration of the parameter space. An excellent discussion of Hamiltonian Monte Carlo is here and the paper extending Hamiltonian Monte Carlo to the No U-Turn Sampler is here. These don't make for light reading and are not required for this course.

The idea: rather than blindly stumbling around the posterior, use the posterior gradiant to skate around the gradiant contour. As you skate closer to a drop-off (gradiant is steep and probability is lower), potential energy decreases and kinetic energy increases (since energy is always conserved). When this happens the skater is turned back uphill and pushed from the precipice and skates on along a posterior likelihood contour. The No U-Turn sampler keeps skating until the skater tries to turn back towards the original point.

One trick to applying this method is to ensure that proposals lead to accepted samples that satisfy the Ergodic Theorem. The papers above go into that topic. In short, they mostly do.

Using pymc3, we can implement NUTS as [windows and mac/linux users can execute this code]:

This same model can be run using glm (general linear model):

So per unit work (per second of computer time), we explore the posterior much more efficiently using NUTS. Here's the video:

A caveat to NUTS as implemented in PyMC3:

NUTS uses a package called Theano (totally in the background) for compiling your model into faster code. While that is nice, sometimes the error messages you might get are less than illuminating.

The following shows the performance of NUTS (versus MH Random Walk) for a variety of interesting posterior shapes:

EMCEE (Affine Invariant Ensemble Sampler Monte Carlo Markov Chain)¶

Also known as the hammer.

EMCEE, along with NUTs, demonstrate the forefront of current methodologies for constructing MCMC's. EMCEE takes a different approach than the NUTs, which is inspired by the movement of molecules in three dimensions. EMCEE constructs the chain by simultaneously unleashing a number of walkers that share information about the posterior in a way that preserves Ergodicity and the properties we need to use our chain for inference. Unlike all nearly every other sampling method, it can do this using parallel computation making it possible to do MCMC inference on large-scale clusters (supercomputers) or machines having many processor cores.

Earlier in the class, I demonstrated that the njobs option in pymc3 could be used to parallelize pymc3 but it needs to be reiterated that this is constructing 3 independent chains and combining them can be done once chain convergence is assessed and validated. In contrast, EMCEE produces a chain from the get-go that can combine work done across a large number of processor cores. This makes EMCEE a very strong option for computationally expensive problems. At present, this is not implemented in pymc3, unfortunately, but as of March 2018, it is unavailable but work is ongoing.

The mathematics behind EMCEE can be found in the paper here and it has been used extensively in astro-physics. Rather than spend a lot of time with the details, I'll give you an idea of how the sampler works. It should be noted that like NUTs, we are talking about how proposals are determined- once we have a proposal, the accept/reject criteria follows Metropolis-Hastings.

Idea: For a $K$ parameter problem, let there be $W$ walkers. You can think of each walker as being an MCMC chain. Proposals for a walker are found by randomly choosing another walker location and leap-frogging towards it. The length of the jump is determined by the relative likelihood of the two proposal values. A lower likelihood and the walker will leap between its current position and the sampled walker. If higher, it can leap over the sampled walkers position.

Some intuition: In the walking dead, walkers are constantly roaming the landscape looking for victims. When a walker runs into a victim, the other walkers tend to go towards it and wander in that area in a frenzy. If we think of victims as high posterior areas, then the walkers will bounce around good areas almost in a swarm.

Implementing EMCEE in code¶

Unfortunately, EMCEE is not implemented in pymc3, so we have to install it in our Anaconda environment:

conda install emcee

Once installed, the process for setting up the problem is similar to pymc3 in abstract, but differs quite a bit in construction

The estimator returns a Markov Chain (actually the Ensemble Sampler returns 16 chains (number of walkers)). As an aside, this calculates the likelihood function 2,000 x 16 = 32k times. This may be a little overkill for an OLS model.....

The following graphs one of our 16 walkers:

See discussion here: https://dfm.io/emcee/current/user/line/ Notice, using full set of walkers gives us a chain that at least visually does not have the periodicity problems (although there are formal tests), so we can proceed.

Using the same information, the following plots the likelihood (posterior) surface in two dimensions. Each point in the diagrams is an actual parameter value that was included in our chain.