Bayesian Linear Regression.

bayes_linear_regression.py

import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Create some data
n = 100
_b0 = 2
_b1 = 6
x = np.linspace(0, 1, n)
y = _b1*x + _b0 + np.random.randn(n)

plt.plot(x, _b1*x + _b0, label='Truth Line')
plt.scatter(x, y, s=30, label='Data')
plt.legend(loc='best')
plt.savefig('linear_data.png')

niter = 1000
with pm.Model() as model:
    # Define model
    b0 = pm.Normal('b0', mu=0, sd=20)
    b1 = pm.Normal('b1', mu=0, sd=20)
    sigma = pm.Uniform('sigma', lower=0, upper=20)
    mu = b1*x + b0
    likelihood = pm.Normal('y', mu=mu, sd=sigma, observed=y)

    # perform inference
    map_estimate = pm.find_MAP()
    step = pm.NUTS()
    trace = pm.sample(niter, step, map_estimate, progressbar=True)

pm.summary(trace)
df = pm.trace_to_dataframe(trace)

# Create MAP model
y_hat = map_estimate['b1']*x + map_estimate['b0']

# Compute 95% (2*sigma) credible intervals from samples
yfit = df['b0'][:, None] + df['b1'][:, None] * x
sig = 2 * yfit.std(axis=0)

# Plot MAP and 2*sigma range
f, axarr = plt.subplots(2)
axarr[0].plot(x, _b1*x + _b0, label='Truth Line')
axarr[0].fill_between(x, y_hat - sig, y_hat + sig, color='lightgray', alpha=0.8)
axarr[0].scatter(x, y, s=30, label='Data')
axarr[0].plot(x, y_hat, label='MAP estimate')
axarr[0].legend(loc='best')

# Plot MAP and some sample lines
axarr[1].plot(x, _b1*x + _b0, label='Truth Line')
for i in range(100):
    axarr[1].plot(x, df.loc[i]['b1']*x + df.loc[i]['b0'], alpha=0.2, c='red')
axarr[1].scatter(x, y, s=30, label='Data')
axarr[1].plot(x, y_hat, label='MAP estimate')
axarr[1].legend(loc='best')

plt.savefig('linear_fit.png')