# Posterior predictive distributions

When you need the full distribution of new observations (and not just a point estimate or a confidence interval) [posterior_predict()](../reference/GAM.md#whittaker.GAM.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()](../reference/GAM.md#whittaker.GAM.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()](../reference/GAM.md#whittaker.GAM.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](../reference/PosteriorPredictResult.md#whittaker.PosteriorPredictResult) object that stores the full `(n, n_draws)` sample matrix and provides convenience methods for common summaries.


``` python
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.


``` python
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()](../reference/PosteriorPredictResult.md#whittaker.PosteriorPredictResult.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.


``` python
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:


``` python
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()](../reference/PosteriorPredictResult.md#whittaker.PosteriorPredictResult.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.


``` python
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.


``` python
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,
)
```


<style>
  #altair-viz-dcb330b57323427d855fd5948d0e676f.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-dcb330b57323427d855fd5948d0e676f.vega-embed details,
  #altair-viz-dcb330b57323427d855fd5948d0e676f.vega-embed details summary {
    position: relative;
  }
</style>


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](ppc.md) 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()](../reference/GAM.md#whittaker.GAM.posterior_predict)** | **Yes** | **[PosteriorPredictResult](../reference/PosteriorPredictResult.md#whittaker.PosteriorPredictResult) with convenience methods** |

[posterior_predict()](../reference/GAM.md#whittaker.GAM.posterior_predict) is equivalent to `simulate(unconditional=True)` but returns a richer result object with [mean()](../reference/PosteriorPredictResult.md#whittaker.PosteriorPredictResult.mean), [std()](../reference/PosteriorPredictResult.md#whittaker.PosteriorPredictResult.std), [quantile()](../reference/PosteriorPredictResult.md#whittaker.PosteriorPredictResult.quantile), and [interval()](../reference/PosteriorPredictResult.md#whittaker.PosteriorPredictResult.interval) methods. Use [posterior_predict()](../reference/GAM.md#whittaker.GAM.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:


``` python
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,
)
```


<style>
  #altair-viz-a10e99f3b5c44a76a5266ab3d81048fa.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-a10e99f3b5c44a76a5266ab3d81048fa.vega-embed details,
  #altair-viz-a10e99f3b5c44a76a5266ab3d81048fa.vega-embed details summary {
    position: relative;
  }
</style>


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()](../reference/GAM.md#whittaker.GAM.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.


``` python
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:


``` python
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,
)
```


<style>
  #altair-viz-7f85892698aa4bcb9cf6429e3eb41d2f.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-7f85892698aa4bcb9cf6429e3eb41d2f.vega-embed details,
  #altair-viz-7f85892698aa4bcb9cf6429e3eb41d2f.vega-embed details summary {
    position: relative;
  }
</style>


This asymmetry is a natural consequence of the Poisson distribution and is captured automatically because [posterior_predict()](../reference/GAM.md#whittaker.GAM.posterior_predict) draws from the actual response distribution rather than relying on a symmetric normal approximation.


# Works with all inference methods

[posterior_predict()](../reference/GAM.md#whittaker.GAM.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.


``` python
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.


``` python
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:


``` python
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](ppc.md).


# Where to go next

- **[Posterior predictive checks](ppc.md)**: automated PPC with built-in test statistics that compare replicated datasets against the observations.
- **[Prediction and inference](prediction.md)**: confidence intervals, simultaneous bands, and term-level predictions for when you need uncertainty in the mean rather than in new observations.
- **[Variational inference](variational-inference.md)**: fitting Bayesian GAMs whose posterior feeds into [posterior_predict()](../reference/GAM.md#whittaker.GAM.posterior_predict).
- **[MCMC sampling](mcmc.md)**: exact posterior inference for the most accurate predictive distributions.
- **[Model comparison with LOO](loo.md)**: comparing models by out-of-sample predictive accuracy using PSIS-LOO.