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 8.7  Bayesian random intercept binary logistic model in Python using pymc3

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

import numpy as np
import pymc3 as pm

from scipy.stats import norm, uniform, bernoulli

â€‹

# Data
np.random.seed(13531)                                                 # set seed to replicate example
N = 4000                                                                       # number of obs in model
NGroups = 20

â€‹

x1 = uniform.rvs(size=N)
x2 = uniform.rvs(size=N)
Groups = np.array([200 * [i] for i in range(20)]).flatten()

â€‹

a = norm.rvs(loc=0, scale=0.5, size=NGroups)
eta = 1 + 0.2 * x1 - 0.75 * x2 + a[list(Groups)]
mu = 1.0/(1.0 + np.exp(-eta))
y = bernoulli.rvs(mu, size=N)

â€‹

with pm.Model() as model:
# Define priors
sigma = pm.Uniform('sigma', 0, 100)
sigma_a = pm.Uniform('sigma_a', 0, 10)
beta1 = pm.Normal('beta1', 0, sd=100)
beta2 = pm.Normal('beta2', 0, sd=100)
beta3 = pm.Normal('beta3', 0, sd=100)

# priors for random intercept (RI) parameters
a_param = pm.Normal('a_param',
np.repeat(0, NGroups),                           # mean
sd=np.repeat(sigma_a, NGroups),          # standard deviation
shape=NGroups)                                     # number of RI parameters

â€‹

eta = beta1 + beta2 * x1 + beta3  *  x2 + a_param[Groups]

# Define likelihood
y = pm.Normal('y', mu=1.0/(1.0 + np.exp(-eta)), sd=sigma, observed=y)

# Fit
start = pm.find_MAP()                                                    # Find starting value by optimization
step = pm.NUTS(scaling=start)                                       # Initiate sampling
trace = pm.sample(7000, step, start=start, progressbar=False)

# Print summary to screen
pm.summary(trace)

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

Output on screen:

â€‹

sigma_interval:

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

-5.397           0.011            0.001            [-5.419, -5.375]

â€‹

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|
-5.418         -5.405         -5.398         -5.389         -5.374

sigma_a_interval:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------
-2.842           0.205            0.012            [-3.293, -2.461]

â€‹

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|
-3.243         -2.976         -2.845         -2.718         -2.399

beta1:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------
1.140            0.160            0.004            [0.824, 1.450]

â€‹

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|
0.831          1.032          1.137          1.245          1.461

beta2:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------
0.223            0.121            0.001            [-0.015, 0.459]

â€‹

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|
-0.014         0.143          0.224          0.305          0.461

beta3:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------
-0.861           0.126            0.001            [-1.125, -0.629]

â€‹

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|
-1.112         -0.945         -0.860         -0.776         -0.612

a_param:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------
-0.343           0.183            0.004            [-0.690, 0.024]
-0.295           0.183            0.004            [-0.671, 0.042]
-0.192           0.187            0.004            [-0.552, 0.181]
-0.341           0.184            0.004            [-0.710, 0.014]
0.103            0.193            0.004            [-0.262, 0.489]
-0.410           0.181            0.004            [-0.766, -0.058]
0.004            0.193            0.004            [-0.382, 0.376]
-0.420           0.181            0.004            [-0.784, -0.066]
-1.063           0.180            0.004            [-1.432, -0.724]
-0.724           0.182            0.004            [-1.074, -0.359]
-0.047           0.187            0.004            [-0.379, 0.364]
0.345            0.205            0.004            [-0.053, 0.746]
0.395            0.210            0.004            [-0.006, 0.817]
0.890            0.250            0.005            [0.394, 1.385]
0.352            0.207            0.004            [-0.044, 0.754]
0.020            0.195            0.004            [-0.357, 0.409]
0.109            0.198            0.004            [-0.266, 0.510]
0.632            0.224            0.004            [0.201, 1.088]
0.859            0.249            0.005            [0.398, 1.377]
0.197            0.199            0.004            [-0.172, 0.601]

â€‹

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|
-0.698         -0.464         -0.343         -0.220         0.018
-0.656         -0.418         -0.292         -0.170         0.062
-0.560         -0.320         -0.192         -0.065         0.175
-0.708         -0.459         -0.340         -0.220         0.018
-0.269         -0.026         0.100          0.229          0.486
-0.764         -0.531         -0.410         -0.289         -0.055
-0.367         -0.126         0.003          0.131          0.393
-0.787         -0.538         -0.417         -0.299         -0.067
-1.426         -1.181         -1.062         -0.941         -0.716
-1.083         -0.847         -0.721         -0.602         -0.362
-0.414         -0.172         -0.052         0.078          0.331
-0.042         0.205          0.338          0.477          0.765
-0.008         0.252          0.390          0.530          0.815
0.419          0.720          0.881          1.052          1.418
-0.034         0.211          0.348          0.487          0.770
-0.356         -0.109         0.018          0.147          0.410
-0.269         -0.022         0.104          0.239          0.509
0.206          0.480          0.625          0.778          1.100
0.401          0.688          0.845          1.017          1.385
-0.191         0.062          0.198          0.331          0.588

sigma:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------
0.451            0.005            0.000            [0.441, 0.461]

â€‹

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

0.441          0.447          0.451          0.455          0.461

sigma_a:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------
0.561            0.112            0.007            [0.353, 0.777]

â€‹

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|
0.376          0.485          0.549          0.619          0.833

bottom of page