Exercises

import shutil
import subprocess

if shutil.which("uv") is None:
    raise RuntimeError(
        "uv is not installed.\n"
        "Install it from https://docs.astral.sh/uv/getting-started/installation/\n"
        "Then restart the kernel and run this cell again."
    )

subprocess.run(["uv", "sync"], check=True)

print("Environment created successfully.")
print("Now change the kernel to uv: stancon2026-workflow-tools")

import arviz as az

Part 1: Posterior draws objects

For this set of exercises we will use draws from the classic eight schools model. These draws are included in both posterior and ArviZ.

dt = az.load_arviz_data('non_centered_eight')
eight_schools_draws = dt.posterior

Understanding draws objects

Inspect the posterior object. How many chains, iterations and variables does it contain?

# Your code here

Subsetting draws

Extract only the first chain.

You can use the sel() method, specifying the chain argument.

# Your code here

Extract only the first 10 iterations, from all the chains.

You can use the sel() method, specifying the draw argument and using the slice() function.

# Your code here

Thinning draws

Thin the draws so that only half of the draws are included. Then try automatic thinning.

You can use the function az.thin(), optionally specifying the factor argument.

# Your code here

Summarising draws

Extract the variable mu and summarise it via mean and sd.

You can use the functions az.summary(), specifying the var_names argument and the kind argument.

# Your code here

Creating new variables

Create a new variable that is the difference between school 1 and school 2 means. Call it diff_1_2. Do the same with school 3 and school 4. Then summarise these new variables with median, 0.3 and 0.7 quantiles.

# Your code here

Marginal posteriors

Plot the marginal posteriors for the school means.

# Your code here

Pairs plot

Plot the variables mu and tau in an pairs plot.

# Your code here

Part 2: Convergence diagnostics and uncertainty

R-hat

Calculate the R-hat for all the variables in the model. Which variables have high R-hat (> 1.01)?

You can use the function az.rhat()

# Your code here

Effective sample size (ESS)

Calculate the bulk and tail ESS for all the variables in the model.

# Your code here

Monte Carlo standard error

Calculate the mean of each variable, and also the Monte Carlo standard error of the mean.

# Your code here

Then do the same for the 0.05 and 0.95 quantiles. Think about how the Monte Carlo standard error might influence how you report the quantiles.

# Your code here

Pareto diagnostics

Calculate the minimum sample size for stable estimates for each variable in the model. Which has the highest minimum sample size?

# Your code here

Part 3: Model evaluation and critique

We can generate prior predictions from the eight schools model, using the following function.

import numpy as np
import pymc as pm
import pymc.dims as pmd

def make_eight_schools_model(mu_prior_sd=1, tau_prior_sd=1, observe=False):
  coords = {"school": eight_schools_draws.coords["school"].values}
  sigma_obs=xr.DataArray(np.array([15, 10, 16, 11, 9, 11, 10, 18]), coords=coords, dims=("school"))
  
  with pm.Model(coords=coords) as model:
    mu = pmd.Normal("mu", mu=0, sigma=mu_prior_sd)
    tau = pmd.HalfNormal("tau", sigma=tau_prior_sd)
    theta = pmd.Normal("theta", mu=mu, sigma=tau, dims=("school"))
    yrep = pmd.Normal("yrep", mu=theta, sigma=sigma_obs, dims=("school"), observed= dt["observed_data"]["obs"] if observe else None)
    
  return model

def eight_schools_prior(ndraws=1000, mu_prior_sd=1, tau_prior_sd=1):
  model = make_eight_schools_model(mu_prior_sd=mu_prior_sd, tau_prior_sd=tau_prior_sd)
  with model:
    prior_predictive_draws = pm.sample_prior_predictive(draws=ndraws)
  return prior_predictive_draws

Generate 1000 prior predictive draws, and plot the distributions for each school. Try with different mu_prior_sd and tau_prior_sd values (e.g. 1, 10, 100).

prior_predictive_draws = eight_schools_prior(
  ndraws=1000,
  mu_prior_sd=1,
  tau_prior_sd=1
)

# Your code here

Posterior predictive checks

We can create posterior predictive draws from our posterior draws and plot against our actual observations.

with make_eight_schools_model() as model:
    posterior_predictive_draws = pm.sample_posterior_predictive(eight_schools_draws, var_names=["yrep"])

posterior_predictive_draws

Plot the posterior predictions on top of the actual observations. Then use the PIT-ECDF plot.

You can use the az.plot_ppc_pit() function.

y = np.array([28,  8, -3,  7, -1,  1, 18, 12])
# Your code here

Sensitivity checks

Check for prior and likelihood sensitivity in the eight schools model. First check by power-scaling all priors jointly, then select only the “mu” and only the “tau” prior separately.

eight_schools_draws
with make_eight_schools_model(observe=True) as model:
  pm.stats.compute_log_prior(
    dt,
    var_names=["mu", "tau", "theta"],
  #  extend_inferencedata=False
    )
# Your code here

Next plot sensitivity using as density plots. Plot only the mu and tau variables.

# Your code here

Part 4: Bringing it all together

We have provided four sets of posterior draws from hierarchical models of observed migratory bird counts recorded between 2000 and 2020 at the Hanko Bird Observatory (Halias).

For species \(j\),

\(y \sim \mathrm{Poisson}(\lambda_j)\)

\(y \sim \mathrm{NegativeBinomial}(\lambda_j, \phi).\)

The species-specific abundances are linked through a hierarchical prior,

\(\log(\lambda_j) \sim \mathrm{Normal}(\mu, \sigma)\)

where \(\mu\) represents the average abundance across species and \(\sigma\) controls the amount of pooling between species.

Your task is to explore the posterior draws and diagnostic outputs for the four fitted models.

As you work through the diagnostics, try to identify which model corresponds to each of the following situations:

Convergence issues caused by poor chain mixing (for example, insufficient warmup).
Inadequate fit to the data caused by the choice of observation model, shown by posterior predictive checks.
Issues caused by priors in conflict with the likelihood.
No major issues, although there is still room for model improvement.

Use posterior summaries, convergence diagnostics, posterior predictive checks, and sensitivity analyses to guide your investigation.

You can also look at the Stan model (birds_per_year.stan) and consider how you might improve it.

# Your code here

Other Formats