# Credible intervals with Capybaras, Python and Stan (cmdstanpy interface)

*R and Shiny Training: If you find this blog to be interesting, please note that I offer personalized and group-based training sessions that may be reserved through Buy me a Coffee. Additionally, I provide training services in the Spanish language and am available to discuss means by which I may contribute to your Shiny project.*

# 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”.

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

```
import pandas as pd
= pd.DataFrame({
baby_capybaras "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 {
1]);
n1 ~ binomial(N1, p[2]);
n2 ~ binomial(N2, p[ }
```

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()
= CmdStanModel(stan_file='baby-capybaras.stan')
mod_cmdstanpy
= {
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
= mod_cmdstanpy.sample(
fit_cmdstanpy = baby_capybaras_list,
data = 42,
seed = 4,
chains = 500,
refresh = 1000,
iter_warmup = 10000
iter_sampling
)
fit_cmdstanpy.summary()
5% 50% 95% N_Eff \
Mean MCSE StdDev -63.578200 0.007388 1.013660 -65.602100 -63.267200 -62.610000 18824.6
lp__ 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
p[
/s R_hat
N_Eff12858.3 1.000140
lp__ 1] 25065.3 0.999977
p[2] 27026.2 0.999925 p[
```

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%
= 0.05
alpha = [0 + alpha / 2, 1 - alpha / 2]
percentiles
'p[1]'].quantile([percentiles[0], percentiles[1]])
fit_cmdstanpy.draws_pd()[
0.025 0.662371
0.975 0.885329
1], dtype: float64 Name: p[
```

These results are equivalent to the ones obtained with `cmdstanr`

.