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.18 Logistic model using pymc3

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

import numpy as np
from scipy.stats import bernoulli, uniform, binom
import pymc3 as pm
import pylab as plt
import pandas

â€‹

def invlogit(x):
""" Inverse logit function.

input: scalar
output: scalar
"""

return 1.0 / (1 + np.exp(-x))

â€‹

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

â€‹

x1 = binom.rvs(1, 0.6, size=nobs)
x2 = uniform.rvs(size=nobs)

beta0 = 2.0
beta1 = 0.75
beta2 = -5.0

xb = beta0 + beta1 * x1 + beta2 * x2
exb = 1.0/(1 + np.exp(-xb))                                               # logit link function

by = binom.rvs(1, exb, size=nobs)

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

# Fit
niter = 5000                                                                       # 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
p = invlogit(beta0 + beta1 * x1 + beta2 * x2)
y_obs = pm.Binomial('y_obs', n=np.ones(nobs), p=p, observed=by)

â€‹

# inference
start = pm.find_MAP()
step = pm.NUTS()
trace = pm.sample(niter, step, start, 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.925            0.082            0.002            [1.777, 2.099]

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

1.765          1.870          1.926          1.980          2.090

beta1:

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

0.754            0.072            0.001            [0.614, 0.890]

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

0.615          0.704          0.754          0.803          0.893

beta2:

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

-4.892           0.143            0.004            [-5.158, -4.594]

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

-5.179         -4.988         -4.894         -4.795         -4.613

bottom of page