In [3]:

    

# http://pymc-devs.github.io/pymc3/getting_started/#a-motivating-example-linear-regression
import numpy as np
import pylab as pyl

%matplotlib inline


    

In [2]:

    

# initialize random number generator
np.random.seed(123)


    

In [10]:

    

# true parameter values
alpha, sigma = 1, 1
beta = [1, 2.5]

# size of dataset
size = 100

# predictor variable 
X1 = np.linspace(0, 1, size)
X2 = np.linspace(0, 0.2, size)

# simulate outcome variable
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma


    

In [13]:

    

fig, axes = pyl.subplots(1, 2, sharex=True, figsize=(10,4))

axes[0].scatter(X1, Y)
axes[1].scatter(X2, Y)

axes[0].set_ylabel('Y')
axes[0].set_xlabel('X1')
axes[1].set_xlabel('X2')

pyl.show()


    

In [32]:

    

from pymc3 import Model, Normal, HalfNormal 
from pymc3 import find_MAP, NUTS, sample, traceplot, summary
from scipy import optimize


    

In [16]:

    

basic_model = Model()

with basic_model:
    
    # priors for unknown parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)
    
    # expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2
    
    # likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)


    

In [20]:

    

map_estimate = find_MAP(model=basic_model) #BFGS

print(map_estimate)


    

    




{'alpha': array(0.940580333930751), 'beta': array([ 1.23897098,  0.2477942 ]), 'sigma_log': array(-0.05643453185271788)}

In [22]:

    

map_estimate = find_MAP(model=basic_model, fmin=optimize.fmin_powell)

print(map_estimate)


    

    




{'alpha': array(0.9416977243355756), 'beta': array([ 1.28287086,  0.02152522]), 'sigma_log': array(-0.05679975327072704)}

In [24]:

    

with basic_model:
    
    # obtain starting values via MAP
    start = find_MAP(fmin=optimize.fmin_powell)
    
    # instantiate sampler
    step = NUTS(scaling=start)
    
    # draw 500 posterior samples
    trace = sample(500, step, start=start)


    

    




 [-----------------100%-----------------] 500 of 500 complete in 3.7 sec





    




/usr/local/lib/python2.7/site-packages/theano/gof/cmodule.py:293: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility
  rval = __import__(module_name, {}, {}, [module_name])

In [29]:

    

trace['alpha'][-5:]


    

    
Out[29]:





array([ 0.91830579,  0.97286547,  0.99565443,  0.90862384,  0.82481802])

In [31]:

    

traceplot(trace)
pyl.show()


    

In [33]:

    

summary(trace)


    

    




alpha:
 
  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.956            0.197            0.012            [0.565, 1.328]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.547          0.831          0.955          1.087          1.327


beta:
 
  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.326            1.777            0.152            [-2.521, 4.463]
  -0.286           8.652            0.746            [-16.533, 17.176]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  -2.403         0.128          1.335          2.582          5.029
  -17.822        -5.840         -1.061         5.950          16.252


sigma_log:
 
  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  -0.029           0.080            0.004            [-0.170, 0.145]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  -0.178         -0.082         -0.031         0.014          0.140


sigma:
 
  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.975            0.080            0.004            [0.839, 1.151]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.837          0.922          0.969          1.014          1.150

In [ ]:

    

n = 400