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