top of page

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 5.23 Probit model in Python

===============================================

import numpy as np
from scipy.stats import norm, uniform, bernoulli
import pymc3 as pm
import pylab as plt
import pandas
import theano.tensor as tsr

â€‹

def probit_phi(x):
"""Probit transformation."""
mu = 0
sd = 1
return 0.5 * (1 + tsr.erf((x - mu) / (sd * tsr.sqrt(2))))

â€‹

# Data
np.random.seed(135)                                       # set seed to replicate example
nobs = 5000                                                     # number of obs in model

x1 = uniform.rvs(size=nobs)
x2 = 2 * uniform.rvs(size=nobs)

beta0 = 2.0                                                      # coefficients for linear predictor
beta1 = 0.75
beta2 = -1.25

xb = beta0 + beta1 * x1 + beta2 * x2              # linear predictor
exb = 1 - norm.sf(xb)                                      # inverse probit link
py = bernoulli.rvs(exb)

df = pandas.DataFrame({'x1': x1, 'x2':x2, 'by': py})       # re-write data

# Fit
niter = 10000                                              # parameters for MCMC

with pm.Model() as model_glm:
# define priors
beta0 = pm.Flat('beta0')
beta1 = pm.Flat('beta1')
beta2 = pm.Flat('beta2')

# define likelihood
theta_p = beta0 + beta1*x1 + beta2 * x2
theta = probit_phi(theta_p)
y_obs = pm.Bernoulli('y_obs', p=theta, observed=py)

# inference
start = pm.find_MAP()                         # find starting value by optimization
step = pm.NUTS()
trace = pm.sample(niter, step, start, random_seed=135, progressbar=True)

â€‹

# Print summary to screen
pm.summary(trace)

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

===============================================

Output on screen:

â€‹

beta0:

Mean             SD               MC Error         95% HPD interval
------------------------------------------------------------------------------

1.980            0.070            0.001            [1.841, 2.111]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

1.846          1.932          1.980          2.029          2.118

beta1:

Mean             SD               MC Error         95% HPD interval
------------------------------------------------------------------------------

0.757            0.081            0.001            [0.599, 0.910]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

0.601          0.700          0.758          0.811          0.914

beta2:

Mean             SD               MC Error         95% HPD interval
------------------------------------------------------------------------------

-1.236           0.046            0.001            [-1.326, -1.147]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

-1.327         -1.268         -1.235         -1.205         -1.148

bottom of page