Case study: the Nile flow, Bayesian inference beyond i.i.d.

1. The three ways data stop being i.i.d.

The sampler evaluates a per-row log-density and sums the rows. Under independence that sum is the log-likelihood, whether or not the rows share one distribution; under dependence the same sum is still the exact log-likelihood as long as each row is the conditional density of its observation given the previous ones. That single idea organises the three regimes.

1. Independent, not identically distributed. Each row has its own law through its covariates. Every regression in the earlier case studies is already here; the theory changes name (Lindeberg-Feller and local asymptotic normality instead of the i.i.d. central limit theorem), the code does not change at all.
2. Identically distributed, not independent. A stationary series with Markov dependence. The joint density factorises as \(p(y_{1:n}) = p(y_{1:p}) \prod_{t > p} p(y_t \mid y_{t-1}, \dots, y_{t-p})\), so the engine’s row sum computes the conditional log-likelihood when the rows carry the lagged values. gpum_ts_model() builds that design.
3. Neither independent nor identically distributed. Dependence and covariates at once: the conditional density of each row also moves with exogenous variables. The Nile with the dam indicator is exactly this.

The asymptotics that license the posterior summaries are worth stating once, plainly. In regime 2 and 3 the score of the conditional factorisation is a martingale difference sequence, so the martingale central limit theorem (Hall and Heyde 1980) plays the role the classical CLT plays for independent data; maximum likelihood for ergodic Markov processes behaves as in Billingsley (1961); and the Bernstein-von Mises theorem carries over to ergodic Markov processes (Borwanker, Kallianpur and Prakasa Rao 1971), which is what lets the credible intervals below be read the usual way. Stationarity and ergodicity are model assumptions; no diagnostic below can prove them.

2. Regime 2: a stationary AR(1) for the flow

The AR(1) conditional density of one row given its lag is Gaussian, and the whole model is three parameters: the level mu, the persistence phi, the innovation scale on the log scale ls. The formula is the conditional log-density of the current value y given y_lag1; the normalising constant is included so the model can enter predictive comparisons later.

m_ar1 <- gpum_ts_model(
  ~ -0.5 * ((y - mu - phi * (y_lag1 - mu)) / exp(ls))^2 - ls -
    0.9189385332046727,
  params = c("mu", "phi", "ls"),
  series = list(y = y), order = 1,
  prior = ~ -0.5 * ((mu - 900) / 500)^2 - 0.5 * (phi / 2)^2 -
    0.5 * ((ls - 5) / 3)^2
)
m_ar1
#> <gpum_ts>
#>   series     : y (n = 100)
#>   Markov order: 1 (conditional rows: 99)
#>   parameters : mu, phi, ls
#>   factorisation: conditional on the first `order` observations

The conditional factorisation gives 99 rows from 100 observations: the first year is conditioned on, the standard treatment of the initial state. The fit is a regular gpu_metropolis() call; the gradient path is the right tool since the log-density is smooth.

fit_ar1 <- gpu_metropolis(m_ar1$model, data = m_ar1$data,
                          n_iter = 30000, n_chains = 8,
                          method = "mala", warmup = "auto",
                          seed = 1, backend = "cpu")
gpum_diagnose(fit_ar1, plot = FALSE)
#> <gpum_diagnose: Converged>
#>   method mala, backend cpu, chains 8, iterations 19125 (warmup 10875, adaptive)
#> 
#>  parameter     mean       sd     q2.5      q50    q97.5   Rhat        ESS
#>         mu 912.7286 32.94144 846.5802 912.9982 977.2109 1.0001 51652.1266
#>        phi   0.5198  0.09018   0.3447   0.5190   0.6990 1.0001 49324.4030
#>         ls   4.9924  0.07211   4.8556   4.9906   5.1390 1.0001 55923.1730
#>       MCSE
#>  0.1449434
#>  0.0004060
#>  0.0003049
#> 
#>   R-hat below 1.01 and ESS at or above 400 in every parameter.

The classical reference for this exact conditional likelihood is arima(), and the Bernstein-von Mises theorem above says the posterior mean and standard deviation must land on the maximum likelihood estimate and its standard error. They do:

ml <- arima(y, order = c(1, 0, 0), method = "ML")
comp <- rbind(
  posterior = c(mean(fit_ar1$draws[, , "phi"]),
                sd(fit_ar1$draws[, , "phi"]),
                mean(fit_ar1$draws[, , "mu"])),
  arima = c(coef(ml)["ar1"], sqrt(vcov(ml)["ar1", "ar1"]),
            coef(ml)["intercept"])
)
colnames(comp) <- c("phi", "se(phi)", "mu")
round(comp, 4)
#>              phi se(phi)       mu
#> posterior 0.5198  0.0902 912.7286
#> arima     0.5063  0.0867 919.5499

The joint posterior view and the information-bound reference read as in the i.i.d. case studies; the gpum_crlb() documentation states why the observed information of the conditional factorisation is the right curvature here (the martingale central limit theorem stands in for Lindeberg-Feller).

gpum_pairs(fit_ar1)

3. Regime 3: the dam as a covariate, dependence plus heterogeneity

The AR(1) treats the 1898 level shift as persistence it has to absorb. The honest model gives the shift its own parameter: an exogenous indicator dam that is 0 through 1898 and 1 after, entering the conditional mean of each row. The rows are now neither independent (the lag is still there) nor identically distributed (the mean moves with the covariate): the fully non-i.i.d. regime, and still one formula.

dam <- as.numeric(year >= 1899)
m_dam <- gpum_ts_model(
  ~ -0.5 * ((y - mu - b * dam - phi * (y_lag1 - mu - b * dam_lag)) /
              exp(ls))^2 - ls - 0.9189385332046727,
  params = c("mu", "b", "phi", "ls"),
  series = list(y = y), order = 1,
  covariates = list(dam = dam, dam_lag = c(0, dam[-length(dam)])),
  prior = ~ -0.5 * ((mu - 1000) / 500)^2 - 0.5 * (b / 500)^2 -
    0.5 * (phi / 2)^2 - 0.5 * ((ls - 5) / 3)^2
)
fit_dam <- gpu_metropolis(m_dam$model, data = m_dam$data,
                          n_iter = 30000, n_chains = 8,
                          method = "mala", warmup = "auto",
                          seed = 1, backend = "cpu")
gpum_diagnose(fit_dam, plot = FALSE)
#> <gpum_diagnose: Converged>
#>   method mala, backend cpu, chains 8, iterations 19125 (warmup 10875, adaptive)
#> 
#>  parameter      mean       sd       q2.5       q50     q97.5   Rhat        ESS
#>         mu 1096.3060 31.21862 1034.66819 1096.2297 1158.2527 1.0002 51445.3877
#>          b -246.8079 36.36949 -319.22515 -246.7005 -175.4527 1.0001 50200.3218
#>        phi    0.1866  0.10396   -0.01659    0.1860    0.3920 1.0001 52686.5840
#>         ls    4.8517  0.07252    4.71468    4.8499    4.9983 1.0001 43964.2427
#>       MCSE
#>  0.1376387
#>  0.1623244
#>  0.0004529
#>  0.0003459
#> 
#>   R-hat below 1.01 and ESS at or above 400 in every parameter.

Two questions, answered from the draws with the estimation-based tools of the Bayesian layer. Did the dam lower the flow, and by how much? And how much persistence remains once the shift is modelled?

gpum_hypothesis(fit_dam, "b", lower = -Inf, upper = 0)
#> <gpum_hypothesis>
#>   parameter b, interval (-Inf, 0)
#>   P(in interval) = 1.0000
#>   P(>= 0)       = 0.0000
hdi(as.vector(fit_dam$draws[, , "b"]))
#>     lower     upper 
#> -318.9583 -175.2337
round(c(phi_ar1 = mean(fit_ar1$draws[, , "phi"]),
        phi_dam = mean(fit_dam$draws[, , "phi"])), 3)
#> phi_ar1 phi_dam 
#>   0.520   0.187

The persistence drops once the level shift is explained: a large part of what the plain AR(1) read as correlation between consecutive years was the two levels of the series. This is the classic lesson of confounded dependence and heterogeneity, visible in four lines.

4. Model comparison that respects time

WAIC and PSIS-LOO ask a leave-one-out question, and leaving out a past observation that later rows condition on is not a forecasting question: their pointwise exchangeability assumption fails on dependent rows. The tool that respects the arrow of time is leave-future-out cross-validation (Burkner, Gabry and Vehtari 2020): refit on the first t - 1 rows, score the one-step-ahead predictive density of row t, slide forward. gpum_lfo() does this exactly, refitting at every step; the sampler is fast enough that no importance-sampling approximation is needed.

lfo_ar1 <- gpum_lfo(m_ar1, L = 79L, n_iter = 4000, n_chains = 4,
                    seed = 1, backend = "cpu")
lfo_dam <- gpum_lfo(m_dam, L = 79L, n_iter = 4000, n_chains = 4,
                    seed = 1, backend = "cpu")
lfo_ar1
#> <gpum_lfo> exact leave-future-out cross-validation
#>   evaluated points : 20 (training window from 79 rows)
#>   elpd_lfo         : -126.18 (se 1.71)
lfo_dam
#> <gpum_lfo> exact leave-future-out cross-validation
#>   evaluated points : 20 (training window from 79 rows)
#>   elpd_lfo         : -125.11 (se 3.14)

The dam model forecasts the held-out final fifth of the series better, which is the verdict the level shift already announced. Both formulas carry their Gaussian normalising constants, so the elpd values are comparable.

5. Observed against generated, in time

The posterior predictive check for a series compares the observed data with trajectories simulated from the posterior: draw parameters, then simulate the AR recursion forward from the first observation. The distribution overlay is the same discipline as in the i.i.d. case studies, gpum_density_compare() on the pooled values, plus a trajectory band in time.

set.seed(2)
draws <- matrix(fit_dam$draws, ncol = fit_dam$n_params)
colnames(draws) <- m_dam$model$params
sim_one <- function(th) {
  ys <- numeric(length(y)); ys[1] <- y[1]
  for (t in 2:length(y)) {
    mt <- th["mu"] + th["b"] * dam[t]
    mtl <- th["mu"] + th["b"] * dam[t - 1]
    ys[t] <- mt + th["phi"] * (ys[t - 1] - mtl) + rnorm(1, 0, exp(th["ls"]))
  }
  ys
}
sims <- sapply(sample(nrow(draws), 200), function(i) sim_one(draws[i, ]))
gpum_density_compare(y, as.vector(sims),
                     main = "Nile flow: observed vs posterior-simulated")

qs <- apply(sims, 1, quantile, probs = c(0.05, 0.5, 0.95))
plot(year, y, type = "l", lwd = 1.5, ylim = range(qs, y),
     main = "Observed series inside the posterior predictive band",
     xlab = "year", ylab = "flow")
polygon(c(year, rev(year)), c(qs[1, ], rev(qs[3, ])),
        col = grDevices::adjustcolor("steelblue", 0.25), border = NA)
lines(year, qs[2, ], col = "steelblue", lwd = 1.5, lty = 2)
lines(year, y, lwd = 1.5)
abline(v = 1898, lty = 3)
legend("topright", bty = "n", lwd = c(1.5, 1.5), lty = c(1, 2),
       col = c("black", "steelblue"),
       legend = c("observed", "posterior median (90% band)"))

6. Where this stops, honestly

The conditional-factorisation pattern covers Markov dependence of any finite order, with or without exogenous covariates, and everything independent. It does not cover models whose likelihood needs a latent recursion (GARCH variances, state-space filters) or a fully coupled covariance (Gaussian processes, y ~ N(mu, K(theta))), because the per-row DSL has no recursion and no matrix operations. For those, Stan’s transformed parameters, sequential Monte Carlo, or INLA are the right tools, and the roadmap’s matrix-DSL tier is where this package would meet them.

References

• Billingsley, P. (1961). Statistical Inference for Markov Processes. University of Chicago Press.
• Borwanker, J., Kallianpur, G., Prakasa Rao, B. L. S. (1971). The Bernstein-von Mises theorem for Markov processes. Annals of Mathematical Statistics, 42(4), 1241-1253.
• Burkner, P.-C., Gabry, J., Vehtari, A. (2020). Approximate leave-future-out cross-validation for Bayesian time series models. Journal of Statistical Computation and Simulation, 90(14), 2499-2523.
• Hall, P., Heyde, C. C. (1980). Martingale Limit Theory and Its Application. Academic Press.
• van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press.