From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida  

 

you are kindly asked to include the complete citation if you used this material in a publication

Code 3.4 Bayesian normal linear model in Python.
==================================================

import numpy as np
from pymc3 import Model, sample, summary, traceplot
from pymc3.glm import glm
import pylab as plt
import pandas

from scipy.stats import uniform, norm

# Data
np.random.seed(1056)                                                        # set seed to replicate example
nobs = 250                                                                          # number of obs in model 
x1 = uniform.rvs(size=nobs)                                              # random uniform variable

beta0 = 2.0                                                                          # intercept
beta1 = 3.0                                                                          # angular coefficient

xb = beta0 + beta1 * x1                                                      # linear predictor
y = norm.rvs(loc=xb, scale=1.0, size=nobs)                      # create y as adjusted


# Fit
df = pandas.DataFrame({'x1': x1, 'y': y})                          # rewrite data


with Model() as model_glm:
    glm('y ~ x1', df)
    trace = sample(5000)


# Output
summary(trace)


# Show graphical output
traceplot(trace)
plt.show()
==================================================

Python 2.7 - Output on screen:

Intercept:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  2.040            0.133            0.003            [1.792, 2.284]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.794          1.955          2.039          2.128          2.287


x1:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  3.134            0.226            0.005            [2.707, 3.574]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  2.690          2.989          3.134          3.285          3.569


sd_log:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  -0.054           0.045            0.001            [-0.140, 0.037]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  -0.140         -0.085         -0.056         -0.025         0.038


sd:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.948            0.043            0.001            [0.865, 1.033]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.869          0.918          0.946          0.975          1.038