# Compound Steps in Sampling¶

This notebook explains how the compound steps work in pymc3.sample function when sampling multiple random variables. We are going to answer the following questions associated with compound steps:

• How does compound steps work?
• What will happen when PyMC3 assign step methods by default?
• How to specify the step methods? What is the order to apply the step methods at each iteration? Is there a way to specify the order of the step methods?
• What are the issues with mixing discrete and continuous samplers, especially with HMC/NUTS?
In [1]:

import pymc3 as pm
import numpy as np
import theano


## Compound steps¶

When sampling a model with multiple free random variables, the compound steps will be needed in the pm.sample function. When the compound steps are involved, the function takes a list of step to generate a list of methods for different random variables. So for example in the following code,

with pm.Model() as m:
rv1 = ... # random variable 1 (continuous)
rv2 = ... # random variable 2 (continuous)
rv3 = ... # random variable 3 (categorical)
...
step1 = pm.Metropolis([rv1, rv2])
step2 = pm.CategoricalGibbsMetropolis([rv3])
trace = pm.sample(..., step=[step1, step2]...)


the compound step is now contain a list of methods. And at each sampling step, it iterates each method, which takes a point as input, and generates a new point as output. The new point is proposed within each step via a stochastic kernel, and if the proposal was rejected by the Metropolis-Hastings criteria it just outputs the original input point.

## Compound steps by default¶

When we call pm.sample(), PyMC3 assigns the best step method to each of the free random variables. Take the following example:

In [2]:

n_ = theano.shared(np.asarray([10, 15]))
with pm.Model() as m:
p = pm.Beta('p', 1., 1.)
ni = pm.Bernoulli('ni', .5)
k = pm.Binomial('k', p=p, n=n_[ni], observed=4)
trace = pm.sample()

Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [p]
>BinaryGibbsMetropolis: [ni]
100%|██████████| 1000/1000 [00:01<00:00, 657.44it/s]
INFO (theano.gof.compilelock): Waiting for existing lock by process '85686' (I am process '85687')
INFO (theano.gof.compilelock): To manually release the lock, delete /Users/yuanhuang/.theano/compiledir_Darwin-17.7.0-x86_64-i386-64bit-i386-3.6.5-64/lock_dir
INFO (theano.gof.compilelock): Waiting for existing lock by process '85686' (I am process '85688')
INFO (theano.gof.compilelock): To manually release the lock, delete /Users/yuanhuang/.theano/compiledir_Darwin-17.7.0-x86_64-i386-64bit-i386-3.6.5-64/lock_dir
INFO (theano.gof.compilelock): Waiting for existing lock by process '85687' (I am process '85688')
INFO (theano.gof.compilelock): To manually release the lock, delete /Users/yuanhuang/.theano/compiledir_Darwin-17.7.0-x86_64-i386-64bit-i386-3.6.5-64/lock_dir


There are two free parameters in the model we would like to sample from, a continuous variable p_logodds__ and a binary variable ni.

In [3]:

m.free_RVs

Out[3]:

[p_logodds__, ni]


When we call pm.sample(), PyMC3 assigns the best step method to each of them. For example, NUTS was assigned to p_logodds__ and BinaryGibbsMetropolis was assigned to ni.

## Specify compound steps¶

But we can also specify the steps manually:

In [4]:

with m:
step1 = pm.Metropolis([p])
step2 = pm.BinaryGibbsMetropolis([ni])


And now we can pass a point to each step, and see what happens. First, let us generate a test point as the input:

In [5]:

point = m.test_point
point

Out[5]:

{'p_logodds__': array(0.), 'ni': array(0)}


Then pass the point to the first step method pm.Metropolis for random variable p.

In [6]:

point, state = step1.step(point=point)
point, state

Out[6]:

({'p_logodds__': array(-0.8881149), 'ni': array(0)},
[{'tune': True, 'accept': 0.7725535007903375}])


As you can see, the value of ni does not change, but p_logodds__ is updated.

And similarly, you can pass the updated point to step2 and get a sample for ni:

In [7]:

point = step2.step(point=point)
point

Out[7]:

{'p_logodds__': array(-0.8881149), 'ni': array(1)}


Compound step works exactly like this by iterating all the steps within the list. In effect, it is a metropolis hastings within gibbs sampling.

Moreover, pm.CompoundStep is called internally by pm.sample(). We can make them explicit as below:

In [8]:

with m:
comp_step1 = pm.CompoundStep([step1, step2])
comp_step1.methods

Out[8]:

[<pymc3.step_methods.metropolis.Metropolis at 0x10bfadba8>,
<pymc3.step_methods.metropolis.BinaryGibbsMetropolis at 0x10d6442b0>]


## Order of step methods¶

When in the default setting, the parameter update order follows the same order of the random variables, and it is assigned automatically. But if you specify the steps, you can change the order of the methods in the list:

In [9]:

with m:
comp_step2 = pm.CompoundStep([step2, step1])
comp_step2.methods

Out[9]:

[<pymc3.step_methods.metropolis.BinaryGibbsMetropolis at 0x10d6442b0>,

In the sampling process, it always follows the same step order in each sample in the Gibbs-like fashion. More precisely, at each update, it iterates over the list of methods where the accept/reject is based on comparing the acceptance rate with $$p \sim \text{Uniform}(0, 1)$$ (by checking whether $$\log p < \log p_{\text {updated}} - \log p_{\text {current}}$$).