Lasso in PyMC3

lasso.py

from pymc import *
from scipy.stats import norm
import pylab as plt

# Same model as the tutorial    
n = 1000    
x1 = norm.rvs(0, 1, size=n)
x2 = -x1 + norm.rvs(0, 10**-3, size=n)
x3 = norm.rvs(0, 1, size=n)

y = 10 * x1 + 10 * x2 + 0.1 * x3
  
with Model() as model:
    # Laplacian prior only works with Metropolis sampler
    coef1 = Laplace('x1', 0, b=1/sqrt(2))
    coef2 = Laplace('x2', 0, b=1/sqrt(2))
    coef3 = Laplace('x3', 0, b=1/sqrt(2))
    
    # Gaussian prior works with NUTS sampler
    #coef1 = Normal('x1', mu = 0, sd = 1)
    #coef2 = Normal('x2', mu = 0, sd = 1)
    #coef3 = Normal('x3', mu = 0, sd = 1)
    
    likelihood = Normal('y', mu= coef1 * x1 + coef2 * x2 + coef3 * x3, tau = 1, observed=y)
    
    #step = Metropolis() # Works just like PyMC2
    start = find_MAP() 
    step = NUTS(scaling = start)
    trace = sample(20000, step, start = start, progressbar=True) 
    trace = trace[1000::5]
    
plt.figure(figsize=(7, 7))
traceplot(trace)
plt.tight_layout()

autocorrplot(trace)
summary(trace)

notebook.ipynb