Motivation

Matsuura’s excellent book Bayesian Statistical Modeling with Stan, R, and Python uses cmdstanr for all the coding examples. You can access the codes from GitHub.

Because I am reading that book, and in the previous post I explained how to obtain a credible interval with rstan and cmdstanr, I will show how to do the same with cmdstanpy.

Problem

We have a dataset about 100 capybaras (hydrochoerus hydrochaeris). Each capybara has 2 babies, one in each season, with each season being “birth1” and “birth2”. The data is:

In this data, the codification is “female = 0” and “male = 1”.

import pandas as pd

baby_capybaras = pd.DataFrame({
    "birth1": [1,0,0,0,1,1,0,1,0,1,0,0,1,1,0,1,1,0,0,0,1,0,0,0,1,0,0,0,0,1,1,1,
               0,1,0,1,1,1,0,1,0,1,1,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0,1,1,0,1,0,0,
               1,0,0,0,1,0,0,1,1,1,1,0,1,0,1,1,1,1,1,0,0,1,0,1,1,0,1,0,1,1,1,0,
               1,1,1,1],
    "birth2": [0,1,0,1,0,1,1,1,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,0,1,1,1,0,1,1,
               1,0,1,0,0,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,
               0,1,1,0,1,1,1,0,0,0,0,0,0,1,0,0,0,1,1,0,0,1,0,0,1,1,0,0,0,1,1,1,
               0,0,0,0]
})

Which is the probability that a capybara has a male baby if the first baby is a female?

Solution

We can calculate posterior beliefs about the conditional probability by using Stan. Before doing computation, let’s assume the births are dependent and that the prior probability is uniform.

The Stan for code for this problem is:

// Input
data {
  int<lower=0> N1;
  int<lower=0> N2;
  int<lower=0> n1;
  int<lower=0> n2;
}

// Parameters
parameters {
  vector<lower=0, upper=1>[2] p;
}

// Model
model {
  n1 ~ binomial(N1, p[1]);
  n2 ~ binomial(N2, p[2]);
}

From Python we use cmdstanpy, which is similar to cmdstanr for R.

With cmdstanpy we do the following:

# conda activate blog
# pip install cmdstanpy

from cmdstanpy import cmdstan_path, install_cmdstan, CmdStanModel

# run once and only once after installing the package
# install_cmdstan()

mod_cmdstanpy = CmdStanModel(stan_file='baby-capybaras.stan')

baby_capybaras_list = {
    'N1': baby_capybaras['birth1'].value_counts()[0],
    'N2': baby_capybaras['birth1'].value_counts()[1],
    'n1': baby_capybaras['birth2'][baby_capybaras['birth1'] == 0].sum(),
    'n2': baby_capybaras['birth2'][baby_capybaras['birth1'] == 1].sum(),
}

# the seed is for reproducibility
# see https://www.independent.co.uk/life-style/history/42-the-answer-to-life-the-universe-and-everything-2205734.html

fit_cmdstanpy = mod_cmdstanpy.sample(
  data = baby_capybaras_list,
  seed = 42, 
  chains = 4,
  refresh = 500,
  iter_warmup = 1000,
  iter_sampling = 10000
)

fit_cmdstanpy.summary()

           Mean      MCSE    StdDev         5%        50%        95%    N_Eff  \
lp__ -63.578200  0.007388  1.013660 -65.602100 -63.267200 -62.610000  18824.6   
p[1]   0.784325  0.000298  0.057167   0.684343   0.788070   0.872173  36695.5   
p[2]   0.415275  0.000338  0.067206   0.306533   0.414058   0.528851  39566.4   

      N_Eff/s     R_hat  
lp__  12858.3  1.000140  
p[1]  25065.3  0.999977  
p[2]  27026.2  0.999925

This implies that the probability of a male birth if the first birth is a female is 0.78.

The code for the 95% credible interval with cmdstanpy is:

# percentiles 0% + alpha/2 and 100% - alpha/2
# in this case 2.5% and 97.5%
alpha = 0.05
percentiles = [0 + alpha / 2, 1 - alpha / 2]

fit_cmdstanpy.draws_pd()['p[1]'].quantile([percentiles[0], percentiles[1]])

0.025    0.662371
0.975    0.885329
Name: p[1], dtype: float64

These results are equivalent to the ones obtained with cmdstanr.