Sometimes, it is very useful to update a set of parameters together. For example, variables that are highly correlated are often good to update together. In PyMC 3 block updating is simple, as example will demonstrate.

Here we have a LASSO regression model where the two coefficients are strongly correlated. Normally, we would define the coefficient parameters as a single random variable, but here we define them separately to show how to do block updates.

First we generate some fake data.



In [10]:

    
%pylab inline
from matplotlib.pylab import *
from pymc3 import * 
import numpy as np 

d = np.random.normal(size=(3, 30))
d1 = d[0] + 4
d2 = d[1] + 4
yd = .2*d1 +.3*d2 + d[2]









    



Populating the interactive namespace from numpy and matplotlib






    



WARNING: pylab import has clobbered these variables: ['prod', 'zeros_like', 'power', 'sinh', 'det', 'sample', 'minimum', 'tan', 'ones_like', 'sin', 'gradient', 'math', 'cosh', 'clip', 'linalg', 'trace', 'dot', 'random', 'ceil', 'where', 'maximum', 'concatenate', 'cos', 'tanh', 'sqrt', 'floor', 'test', 'fft', 'log', 'info', 'sum', 'flatten', 'exp']
`%matplotlib` prevents importing * from pylab and numpy



In [ ]:

Then define the random variables.



In [11]:

    
lam = 3

with Model() as model:
    s = Exponential('s', 1)
    tau = Uniform('tau', 0, 1000)
    b = lam * tau
    m1 = Laplace('m1', 0, b)
    m2 = Laplace('m2', 0, b)
    
    p = d1*m1 + d2*m2
    
    y = Normal('y', mu=p, sd=s, observed=yd)

For most samplers, including Metropolis and HamiltonianMC, simply pass a list of variables to sample as a block. This works with both scalar and array parameters.



In [12]:

    
import pymc3 
pymc3.variational.advi( model=model, n=100000)









    Out[12]:





({'m1': array(0.48459033591027423),
  'm2': array(0.04708045902248652),
  's_log': array(0.2646004713715775),
  'tau_interval': array(-8.314564627550315)},
 {'m1': array(-3.0429711203909964),
  'm2': array(-3.115115750655253),
  's_log': array(-1.932210086480891),
  'tau_interval': array(0.21826592793700966)})



In [13]:

    
with model: 
    start = find_MAP()
    
    step1 = Metropolis([m1, m2])
    
    step2 = Slice([s, tau])
    
    trace = sample(10000, [step1, step2], start=start)









    



 [-----------------100%-----------------] 10000 of 10000 complete in 29.6 sec



In [14]:

    
traceplot(trace);



In [5]:

    
hexbin(trace[m1],trace[m2], gridsize = 50)









    Out[5]:





<matplotlib.collections.PolyCollection at 0xcc8d64c>