Posterior predictive distributions

When you need the full distribution of new observations (and not just a point estimate or a confidence interval) posterior_predict() returns the complete posterior predictive sample. Each draw accounts for both uncertainty in the model coefficients and observation-level noise from the response distribution, so the result reflects what new data from the same process might actually look like.

This is distinct from a confidence or credible interval for the mean, which only tells you where the average response is likely to be. A posterior predictive interval tells you where an individual new observation is likely to fall, which is almost always wider because it adds the inherent randomness of the response on top of the parameter uncertainty.

When to use posterior_predict

Use posterior_predict() when:

  • you need prediction intervals that include observation noise, not just uncertainty in the mean
  • you want to propagate predictive uncertainty through a downstream calculation (e.g., computing the probability that a new observation exceeds a threshold)
  • you are building custom posterior predictive checks beyond the built-in ppc() statistics
  • you want to visualize the full predictive distribution at specific covariate values

For uncertainty in the mean response only (no observation noise), use simulate(unconditional=False) or predict(interval="credible") instead.

Basic usage

We start by fitting a Gaussian GAM with variational inference and then drawing 2000 posterior predictive samples at 200 evenly spaced new points. The result is a PosteriorPredictResult object that stores the full (n, n_draws) sample matrix and provides convenience methods for common summaries.

import numpy as np
import whittaker as wk

rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2 * np.pi, n)
y = np.sin(x) + rng.normal(0, 0.3, n)

data = {"x": x, "y": y}

model = wk.GAM("y ~ s(x)").fit(data, method="VI")

x_new = np.linspace(0, 2 * np.pi, 200)
pp = model.posterior_predict({"x": x_new}, n_draws=2000, seed=0)
pp
PosteriorPredictResult(n_obs=200, n_draws=2000)

The printed representation shows the number of prediction points and the number of draws. The raw sample matrix is available as pp.samples, a NumPy array of shape (200, 2000) where each column is one complete draw from the posterior predictive distribution.

Summarizing the predictive distribution

Working with 2000 draws per point is useful for downstream calculations, but most of the time you want a compact summary. The result object provides convenience methods for the most common ones.

Mean and standard deviation

The posterior predictive mean is a natural point estimate that averages over both coefficient uncertainty and observation noise. The standard deviation quantifies the total predictive spread at each point.

pp_mean = pp.mean()
pp_std = pp.std()

print(f"Predictive mean shape: {pp_mean.shape}")
print(f"Predictive std shape:  {pp_std.shape}")
print(f"Mean std across points: {pp_std.mean():.4f}")
Predictive mean shape: (200,)
Predictive std shape:  (200,)
Mean std across points: 0.3164

Posterior predictive intervals

interval() returns an equal-tailed interval at the requested coverage level. The default is 95%, meaning 2.5% of the draws fall below the lower bound and 2.5% fall above the upper bound at each prediction point.

lower, upper = pp.interval()
print(f"95% interval width (mean): {np.mean(upper - lower):.4f}")
95% interval width (mean): 1.2353

You can request any coverage level. An 80% interval is narrower because it discards more of the tails:

lower_80, upper_80 = pp.interval(0.80)
print(f"80% interval width (mean): {np.mean(upper_80 - lower_80):.4f}")
80% interval width (mean): 0.8107

Arbitrary quantiles

For more fine-grained summaries, quantile() accepts a scalar or a list of quantile values. A scalar returns one value per prediction point; a list returns a matrix with one row per quantile.

median = pp.quantile(0.5)
print(f"Median shape: {median.shape}")

deciles = pp.quantile([0.1, 0.5, 0.9])
print(f"Deciles shape: {deciles.shape}")
Median shape: (200,)
Deciles shape: (3, 200)

Visualizing predictive uncertainty

Plotting the 95% posterior predictive interval alongside the observed data gives an intuitive picture of where the model expects new observations to fall. The band should contain roughly 95% of the data points if the model is well calibrated.

import altair as alt

obs_data = [{"x": float(x[i]), "y": float(y[i])} for i in range(n)]
points = alt.Chart({"values": obs_data}).mark_circle(
    size=15, opacity=0.3, color="steelblue"
).encode(
    x=alt.X("x:Q", title="x"),
    y=alt.Y("y:Q", title="y"),
)

fit_data = [
    {
        "x": float(x_new[i]),
        "mean": float(pp_mean[i]),
        "lower": float(lower[i]),
        "upper": float(upper[i]),
    }
    for i in range(len(x_new))
]

line = alt.Chart({"values": fit_data}).mark_line(
    color="firebrick", strokeWidth=2
).encode(x="x:Q", y="mean:Q")

band = alt.Chart({"values": fit_data}).mark_area(
    opacity=0.15, color="firebrick"
).encode(x="x:Q", y="lower:Q", y2="upper:Q")

(band + line + points).properties(
    title="Posterior predictive: fitted mean and 95% interval",
    width=500,
    height=300,
)

The band covers approximately 95% of new observations. Points falling outside the band are expected (about 5% of the time for a well-calibrated model). If significantly more or fewer points fall outside, the model may be mis-specified so consider running a formal posterior predictive check to investigate.

Posterior predict vs. other uncertainty tools

Whittaker provides several ways to quantify uncertainty. The key distinction is whether the result includes observation noise (the randomness of individual data points around the mean) or only reflects uncertainty in the estimated mean itself. The table below summarizes all the options:

Method Includes observation noise? Returns
predict(interval="confidence") No Interval for the mean
predict(interval="prediction") Yes (normal approx.) Interval for a new observation
predict(interval="credible") No Bayesian interval for the mean
simulate(unconditional=False) No (n, n_sim) draws of the mean
simulate(unconditional=True) Yes (n, n_sim) draws of new observations
posterior_predict() Yes PosteriorPredictResult with convenience methods

posterior_predict() is equivalent to simulate(unconditional=True) but returns a richer result object with mean(), std(), quantile(), and interval() methods. Use posterior_predict() when you want both the full sample and convenient summaries; use simulate() when you only need the raw matrix.

Threshold exceedance probabilities

One of the most powerful uses of the full posterior predictive sample is estimating the probability that a new observation exceeds (or falls below) a given threshold. This is a calculation that point estimates and intervals cannot provide (you need the entire distribution).

For example, you might ask: “at each value of x, what is the probability that a new observation will be above 0.5?” With the sample matrix in hand, this can be determined as such:

threshold = 0.5
prob_above = np.mean(pp.samples > threshold, axis=1)

thresh_data = [
    {"x": float(x_new[i]), "prob": float(prob_above[i])}
    for i in range(len(x_new))
]

alt.Chart({"values": thresh_data}).mark_line(
    color="steelblue", strokeWidth=2
).encode(
    x=alt.X("x:Q", title="x"),
    y=alt.Y("prob:Q", title=f"P(y_new > {threshold})", scale=alt.Scale(domain=[0, 1])),
).properties(
    title=f"Probability that a new observation exceeds {threshold}",
    width=500,
    height=250,
)

The curve tracks the sine wave: the probability is highest where the true function peaks above the threshold and drops to near zero at the troughs. You can replace 0.5 with any threshold relevant to your application (e.g., a regulatory limit, a clinical cutoff, a business target, etc.).

Poisson example

posterior_predict() respects the family’s observation model. For a Poisson GAM the response distribution is discrete, so every draw is a non-negative integer. This is in contrast to predict(interval="prediction"), which uses a normal approximation and can produce non-integer or negative bounds for count data.

from whittaker.families.poisson import Poisson

rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2 * np.pi, n)
lam = np.exp(1.5 * np.sin(x))
y = rng.poisson(lam).astype(float)

pois_model = wk.GAM("y ~ s(x)", family=Poisson()).fit({"x": x, "y": y}, method="VI")

x_pois = np.linspace(0, 2 * np.pi, 200)
pp_pois = pois_model.posterior_predict({"x": x_pois}, n_draws=2000, seed=0)

print(f"Min draw: {pp_pois.samples.min()}")
print(f"All integer-valued: {np.all(pp_pois.samples == np.round(pp_pois.samples))}")
Min draw: 0.0
All integer-valued: True

The posterior predictive interval for the Poisson model is asymmetric. It’s wider where the predicted rate is high (right-skewed count distribution) and tighter near zero where the distribution is compressed against the lower bound:

lower_pois, upper_pois = pp_pois.interval()
mean_pois = pp_pois.mean()

obs_pois = [{"x": float(x[i]), "y": float(y[i])} for i in range(n)]
points_pois = alt.Chart({"values": obs_pois}).mark_circle(
    size=15, opacity=0.3, color="steelblue"
).encode(x=alt.X("x:Q", title="x"), y=alt.Y("y:Q", title="count"))

fit_pois = [
    {
        "x": float(x_pois[i]),
        "mean": float(mean_pois[i]),
        "lower": float(lower_pois[i]),
        "upper": float(upper_pois[i]),
    }
    for i in range(len(x_pois))
]

line_pois = alt.Chart({"values": fit_pois}).mark_line(
    color="firebrick", strokeWidth=2
).encode(x="x:Q", y="mean:Q")

band_pois = alt.Chart({"values": fit_pois}).mark_area(
    opacity=0.15, color="firebrick"
).encode(x="x:Q", y="lower:Q", y2="upper:Q")

(band_pois + line_pois + points_pois).properties(
    title="Poisson GAM: posterior predictive mean and 95% interval",
    width=500,
    height=300,
)

This asymmetry is a natural consequence of the Poisson distribution and is captured automatically because posterior_predict() draws from the actual response distribution rather than relying on a symmetric normal approximation.

Works with all inference methods

posterior_predict() works with every fitting method. For Bayesian fits (method="VI" or method="MCMC"), draws come from the full posterior. For frequentist fits (method="REML", "GCV", "ML"), the Laplace approximation to the posterior is used instead. This means you can switch inference methods without changing any downstream code that consumes the posterior predictive sample.

model_reml = wk.GAM("y ~ s(x)").fit({"x": x, "y": np.sin(x) + rng.normal(0, 0.3, n)})
pp_reml = model_reml.posterior_predict(n_draws=500, seed=0)
print(f"REML posterior predict: {pp_reml}")
REML posterior predict: PosteriorPredictResult(n_obs=300, n_draws=500)

For Gaussian models with moderate-to-large sample sizes, the Laplace approximation is very accurate, so the posterior predictive samples from a frequentist fit will be nearly indistinguishable from those of a VI or MCMC fit. For non-Gaussian families at small sample sizes, VI or MCMC will generally give better-calibrated predictive distributions.

Training data predictions

When new_data is omitted, predictions are made at the training data points. This is useful for in-sample posterior predictive checks or for comparing the observed response against the model’s predictive distribution at each training observation.

pp_train = model.posterior_predict(n_draws=500, seed=0)
print(f"Training data shape: {pp_train.samples.shape}")
Training data shape: (300, 500)

You can use this to compute, for example, the proportion of training observations that fall within the model’s 95% posterior predictive interval. Here’s a quick calibration check:

lower_train, upper_train = pp_train.interval()
y_train = data["y"]
coverage = np.mean((y_train >= lower_train) & (y_train <= upper_train))
print(f"Empirical coverage: {coverage:.1%}")
Empirical coverage: 96.3%

A well-calibrated model should show coverage close to 95%. Substantially lower coverage suggests the model is overconfident (intervals too narrow), while higher coverage suggests it is overly conservative. For a more thorough assessment, use the built-in posterior predictive checks.

Where to go next