Example of Gamma GLM in PyMC2

Example of Gamma GLM, based on simulated data from Sean Anderson's blog

In [1]:
%matplotlib inline
import numpy as np
import pymc as pm
import matplotlib.pyplot as plt

plt.style.use('fivethirtyeight')

Simulate data under shape-rate configuration of the gamma used in PyMC

In [2]:
N = 100
x = pm.runiform(-1,1,N)
a_true = 0.5
b_true = 1.2
y_true = np.exp(a_true+b_true*x)
shape_true = 10
y = pm.rgamma(alpha=shape_true,beta=shape_true/y_true,size=N)
In [3]:
plt.scatter(x,y)
Out[3]:
<matplotlib.collections.PathCollection at 0x10e8833d0>

Build pymc model

In [4]:
# Priors
a = pm.Normal('a', mu=0.0, tau=0.0001)
b = pm.Normal('b', mu=0.0, tau=0.0001)
shape = pm.Uniform('shape', lower=0., upper=100.)

# Model
mu = pm.Lambda('mu', lambda a=a,b=b,shape=shape: np.exp(a+b*x) )

# Likelihood
Yi = pm.Gamma('Yi', alpha=shape, beta=shape/mu, value=y, observed=True)

# Predictive model fit
xnew = np.linspace(min(x),max(x),N)
Ex_mu = pm.Lambda('Ex_mu', lambda a=a,b=b: np.exp(a+b*xnew))
In [5]:
M = pm.MCMC(locals())
M.sample(10000,9000)
 [-----------------100%-----------------] 10000 of 10000 complete in 2.5 sec
In [6]:
plt.hist(M.a.trace())
plt.plot(a_true,1,marker='.',markersize=30)
Out[6]:
[<matplotlib.lines.Line2D at 0x10e91bd90>]
In [7]:
plt.hist(M.b.trace())
plt.plot(b_true,1,marker='.',markersize=30)
Out[7]:
[<matplotlib.lines.Line2D at 0x10ecbafd0>]
In [8]:
plt.hist(M.shape.trace())
plt.plot(shape_true,1,marker='.',markersize=30)
Out[8]:
[<matplotlib.lines.Line2D at 0x10ef2b410>]
In [14]:
# Plot data and model
plt.scatter(x,y,c="dodgerblue",s=30)
ymu = np.array([np.median(i) for i in M.Ex_mu.trace().T])
yl95 = np.array([np.percentile(i,2.5) for i in M.Ex_mu.trace().T])
yu95 = np.array([np.percentile(i,97.5) for i in M.Ex_mu.trace().T])
plt.plot(xnew,ymu,c="red",linewidth=1.5)
plt.plot(xnew,yl95,c="grey",linewidth=1)
plt.plot(xnew,yu95,c="grey",linewidth=1)
Out[14]:
[<matplotlib.lines.Line2D at 0x10fc27950>]