Example of Gamma GLM, based on simulated data from Sean Anderson's blog
%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
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)
plt.scatter(x,y)
<matplotlib.collections.PathCollection at 0x10e8833d0>
Build pymc model
# 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))
M = pm.MCMC(locals())
M.sample(10000,9000)
[-----------------100%-----------------] 10000 of 10000 complete in 2.5 sec
plt.hist(M.a.trace())
plt.plot(a_true,1,marker='.',markersize=30)
[<matplotlib.lines.Line2D at 0x10e91bd90>]
plt.hist(M.b.trace())
plt.plot(b_true,1,marker='.',markersize=30)
[<matplotlib.lines.Line2D at 0x10ecbafd0>]
plt.hist(M.shape.trace())
plt.plot(shape_true,1,marker='.',markersize=30)
[<matplotlib.lines.Line2D at 0x10ef2b410>]
# 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)
[<matplotlib.lines.Line2D at 0x10fc27950>]