8. Cumulative Example*

import pymc as pm
import numpy as np
import arviz as az
import matplotlib.pyplot as plt
import xarray as xr
import pandas as pd

rng = np.random.default_rng(6420)

8. Cumulative Example*#

Adapted from Unit 9: cumulative2.odc.

Notes#

From the OpenBUGS manual:

cumulative(s1, s2): tail area of distribution of s1 up to the value of s2, s1 must be stochastic, s1 and s2 can be the same

Fitting the model#

As we said on the previous page, we want to evaluate the Empirical CDF at the values of our response variable to check how well our model fits the data. First, we’ll fit our model and add posterior predictive samples to the trace.

# fmt: off
y = np.array(
    [
        0.0, 1.0, 2.0, -1.0, 0.4, -0.5, 0.7, -1.2, 0.1, -0.4, 
        0.2, -0.5, -1.4, 1.8, 0.2, 0.3, -0.6, 1.1, 5.1, -6.3
    ]
)
# fmt: on

with pm.Model() as model:
    # Priors
    theta = pm.Flat("theta")
    ls = pm.Flat("ls")
    s = pm.math.exp(ls)

    # Likelihood
    pm.Normal("y", mu=theta, sigma=s, observed=y)

    # Sampling
    trace = pm.sample(2000, random_seed=rng)
    pm.sample_posterior_predictive(trace, extend_inferencedata=True)
    pm.compute_log_likelihood(trace)

Hide code cell output

 
NUTS[nutpie]: [theta, ls]



Sampling: [y]



















az.summary(trace)










  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      theta
      0.05
      0.491
      -0.72
      0.82
      6476
      5327
      1.00
      0.0061
      0.0048
    

    

      ls
      0.752
      0.168
      0.5
      1
      6379
      5447
      1.00
      0.0021
      0.0016
    

  









Plot posterior predictive distribution vs. observed values#


Now we have posterior samples, the equivalent of the stochastic part of the y[i] node in BUGS. There are different ways to compare these with our ECDF. The easiest is the arviz.plot_ppc function with the kind="cumulative" argument.






az.plot_ppc_dist(trace, kind="ecdf", num_samples=100).show()








../_images/1bb95302faf74927a4a8c2e45d478c863936ab2d3a0b480fe52670824e46dd89.png





Doesn’t look like our posterior does a good job of fitting the observed data. The mean fits okay, but both tails have some issues because of some extreme values near the endpoints. This plot is just comparing individual ECDFs, though. It’s not evaluating the posterior ECDF at the observed \(y\) values.






az.plot_loo_pit(trace, var_names="y").show()








../_images/0ffaefbb43e05ebae624f7eff82aa012a865f1acde8a4c5fdfd2891636745ea7.png









loo_pit = az.loo_pit(trace, var_names="y")["y"].to_numpy()

df = pd.DataFrame({"y": y, "loo_pit": loo_pit})

p = np.clip(loo_pit, 1e-10, 1 - 1e-10)
df["squared_logit"] = np.log(p / (1 - p)) ** 2

df.sort_values("squared_logit", ascending=False).style.format(
    "{:.3f}"
).background_gradient(subset=["loo_pit"], cmap="binary", vmin=0, vmax=1)












  
    

       
      y
      loo_pit
      squared_logit
    

  
  
    

      19
      -6.300
      0.000
      70.317
    

    

      18
      5.100
      0.996
      31.098
    

    

      2
      2.000
      0.826
      2.423
    

    

      13
      1.800
      0.798
      1.883
    

    

      12
      -1.400
      0.234
      1.401
    

    

      7
      -1.200
      0.280
      0.893
    

    

      17
      1.100
      0.701
      0.724
    

    

      3
      -1.000
      0.305
      0.679
    

    

      1
      1.000
      0.670
      0.500
    

    

      16
      -0.600
      0.373
      0.270
    

    

      6
      0.700
      0.616
      0.223
    

    

      11
      -0.500
      0.392
      0.191
    

    

      5
      -0.500
      0.408
      0.140
    

    

      9
      -0.400
      0.410
      0.132
    

    

      4
      0.400
      0.571
      0.081
    

    

      15
      0.300
      0.544
      0.032
    

    

      10
      0.200
      0.532
      0.017
    

    

      14
      0.200
      0.524
      0.009
    

    

      0
      0.000
      0.486
      0.003
    

    

      8
      0.100
      0.497
      0.000
    

  







The last two values (indices 18 and 19) are the outliers. Notice the squared logit transformation is useful because it’s sortable by magnitude.






Another cool Arviz plot#


Arviz has the plot_lm function as well. It’s really nice for visualizing the observed data versus the posterior predictive samples. Let’s try it with the single-predictor model from 11. Brozek index prediction*.






data = pd.read_csv("../data/fat.tsv", sep="\t")

y = data["y"].to_numpy(copy=True)
X = data["X8"].to_numpy(copy=True)

with pm.Model():
    x = pm.Data("X8", X, dims="obs")

    sigma = pm.Exponential("sigma", 0.05)
    beta0 = pm.Normal("intercept", 0, sigma=10)
    beta1 = pm.Normal("slope", 0, sigma=10)

    mu = beta0 + beta1 * x
    likelihood = pm.Normal("y", mu=mu, sigma=sigma, observed=y, dims="obs")

    trace = pm.sample(2000, random_seed=rng)
    pm.sample_posterior_predictive(
        trace, extend_inferencedata=True, random_seed=rng
    )











Hide code cell output





NUTS[nutpie]: [sigma, intercept, slope]







Sampling: [y]

















az.plot_lm(trace, x="X8", y="y").show()











Hide code cell output




../_images/06c1f5562bf2b6a8f6c1addaeb3415475dedd74a40f15c45f6ecc74d91c923c6.png











%load_ext watermark
%watermark -n -u -v -iv -p pytensor









Last updated: Tue, 10 Nov 2026

Python implementation: CPython
Python version       : 3.14.7
IPython version      : 9.16.1

pytensor: 3.3.0

arviz     : 1.3.0
matplotlib: 3.11.1
numpy     : 2.4.6
pandas    : 3.0.5
pymc      : 6.3.2
xarray    : 2026.7.0