Skip to main content

Custom Priors

Refine prior distributions using an iterative inspect-check-refine workflow with prior predictive checks.
Author

Benjamin Vincent

Every Bayesian model encodes assumptions about plausible parameter values before seeing data. pathmc ships with weakly informative defaults that work out of the box, but domain knowledge almost always lets you do better. Tighter, more realistic priors improve sampling efficiency, reduce the influence of prior-data conflict, and — critically for causal inference — produce more credible interventional predictions from do().

This example walks through the iterative prior workflow: inspect the defaults, run a prior predictive check, refine, and verify.

Simulate data with known structure

We generate data from a simple causal chain where a treatment X affects an outcome Y through a mediator M, with known coefficients.

X X M M X->M a = 0.5 Y Y X->Y c = 0.3 M->Y b = 0.8
Figure 1: Causal DAG: X affects Y directly (c = 0.3) and indirectly through M (a = 0.5, b = 0.8).
import numpy as np
import pandas as pd
import pathmc
from pathmc import Prior

rng = np.random.default_rng(42)
n = 500

X = rng.normal(size=n)
M = 0.5 * X + rng.normal(scale=0.3, size=n)
Y = 0.8 * M + 0.3 * X + rng.normal(scale=0.3, size=n)

df = pd.DataFrame({"X": X, "M": M, "Y": Y})

True values: a = 0.5, b = 0.8, c = 0.3, \sigma_M = 0.3, \sigma_Y = 0.3.

Inspect the default priors

When you call pathmc.model() without specifying priors, sensible defaults are used. The .equations() method shows the structural equations alongside the priors.

spec = """
M ~ a*X
Y ~ b*M + c*X
"""

model = pathmc.model(spec, data=df)
model.equations()

\begin{aligned} \beta_{M} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{M} &\sim \text{HalfNormal}(sigma=1) \\ \beta_{Y} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{Y} &\sim \text{HalfNormal}(sigma=1) \\[6pt] \mu_{M} &= \beta_{0,\,M} + a \cdot \mathrm{X} \\ \mathrm{M} &\sim \text{Normal}(\mu_{M},\, \sigma_{M}) \\ \mu_{Y} &= \beta_{0,\,Y} + b \cdot \mathrm{M} + c \cdot \mathrm{X} \\ \mathrm{Y} &\sim \text{Normal}(\mu_{Y},\, \sigma_{Y}) \end{aligned}

The defaults are deliberately vague: Normal(mu=0, sigma=10) for regression coefficients and HalfNormal(sigma=1) for residual standard deviations. These are reasonable starting points when you know nothing about the scale of your data, but they may be wider than necessary.

Prior predictive check: are defaults plausible?

A prior predictive check asks: “If we sampled parameters from our priors and generated fake data, would it look remotely like real-world data?” If the prior implies outcomes spanning thousands of units when the real data lives in [-3, 3], the priors are too vague.

Code
import matplotlib.pyplot as plt
import arviz as az

FIG_WIDTH = 7
FIG_HEIGHT = 3.5

COLOR_DEFAULT = "#1b9e77"
COLOR_CUSTOM = "#d95f02"
COLOR_TRUE = "black"
ppc_default = model.sample_prior_predictive(draws=500, random_seed=42)
Sampling: [M, Y, beta_M, beta_Y, sigma_M, sigma_Y]
Code
y_prior = ppc_default.prior["Y"].values.flatten()

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))
x_kde, y_kde, _ = az.kde(y_prior)
ax.plot(x_kde, y_kde, color=COLOR_DEFAULT, lw=2, label="Prior predictive (default)")
ax.fill_between(x_kde, y_kde, alpha=0.6, color=COLOR_DEFAULT)
ax.axvline(
    df["Y"].mean(),
    color=COLOR_TRUE,
    linestyle="--",
    linewidth=1.5,
    label=f"Observed mean = {df['Y'].mean():.2f}",
)
ax.set_xlabel("Y")
ax.set_ylabel("Density")
ax.legend()
ax.set_title("Prior predictive: default priors")
plt.tight_layout()
plt.show()
Figure 2: Prior predictive distribution of Y under default priors. The wide tails (spanning well beyond ±10) indicate the model considers extreme outcomes plausible before seeing data.

The default priors generate Y values spanning a huge range. For this dataset, where Y lives roughly in [-2, 2], coefficients of \pm 10 are implausible. We can do better.

Refine priors with set_priors()

The set_priors() method lets you override specific priors without re-specifying the model. It merges your changes with the current configuration and recompiles.

Since we know the predictors and outcomes are roughly standard-normal, regression coefficients above \pm 3 in magnitude are implausible. We also know the residual noise should be moderate — not zero, but not 10 either.

model.set_priors({
    "beta_M": Prior("Normal", mu=0, sigma=2),
    "beta_Y": Prior("Normal", mu=0, sigma=2),
    "sigma_M": Prior("HalfNormal", sigma=0.5),
    "sigma_Y": Prior("HalfNormal", sigma=0.5),
})

model.equations()

\begin{aligned} \beta_{M} &\sim \text{Normal}(mu=0,\, sigma=2) \\ \sigma_{M} &\sim \text{HalfNormal}(sigma=0.5) \\ \beta_{Y} &\sim \text{Normal}(mu=0,\, sigma=2) \\ \sigma_{Y} &\sim \text{HalfNormal}(sigma=0.5) \\[6pt] \mu_{M} &= \beta_{0,\,M} + a \cdot \mathrm{X} \\ \mathrm{M} &\sim \text{Normal}(\mu_{M},\, \sigma_{M}) \\ \mu_{Y} &= \beta_{0,\,Y} + b \cdot \mathrm{M} + c \cdot \mathrm{X} \\ \mathrm{Y} &\sim \text{Normal}(\mu_{Y},\, \sigma_{Y}) \end{aligned}

ppc_custom = model.sample_prior_predictive(draws=500, random_seed=42)
Sampling: [M, Y, beta_M, beta_Y, sigma_M, sigma_Y]
Code
y_default = ppc_default.prior["Y"].values.flatten()
y_custom = ppc_custom.prior["Y"].values.flatten()

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))

x_kde, y_kde, _ = az.kde(y_default)
ax.plot(x_kde, y_kde, color=COLOR_DEFAULT, lw=2, label="Default priors")
ax.fill_between(x_kde, y_kde, alpha=0.45, color=COLOR_DEFAULT)
x_kde, y_kde, _ = az.kde(y_custom)
ax.plot(x_kde, y_kde, color=COLOR_CUSTOM, lw=2, label="Custom priors")
ax.fill_between(x_kde, y_kde, alpha=0.55, color=COLOR_CUSTOM)
ax.axvline(
    df["Y"].mean(),
    color=COLOR_TRUE,
    linestyle="--",
    linewidth=1.5,
    label=f"Observed mean = {df['Y'].mean():.2f}",
)
ax.set_xlabel("Y")
ax.set_ylabel("Density")
ax.legend()
ax.set_title("Prior predictive comparison")
plt.tight_layout()
plt.show()
Figure 3: Prior predictive distributions of Y under default (green) and refined (orange) priors. The custom priors concentrate mass in the plausible range without excluding any realistic values.

The custom priors concentrate probability mass in the plausible range. The model still allows coefficients up to \pm 4–6 (the 2\sigma range of Normal(0, 2)), which is generous, but it no longer considers coefficients of 10 or 20 as likely.

Sample and verify

With priors that encode our domain knowledge, we can proceed to posterior inference.

model.fit(random_seed=42)
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta_Y, sigma_M, beta_M, sigma_Y]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 3 seconds.
<xarray.DataTree>
Group: /
├── Group: /posterior
│       Dimensions:       (chain: 4, draw: 1000, Y_predictors: 3, M_predictors: 2,
│                          mu_Y_dim_0: 500, mu_M_dim_0: 500)
│       Coordinates:
│         * chain         (chain) int64 32B 0 1 2 3
│         * draw          (draw) int64 8kB 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
│         * Y_predictors  (Y_predictors) <U9 108B 'Intercept' 'M' 'X'
│         * M_predictors  (M_predictors) <U9 72B 'Intercept' 'X'
│         * mu_Y_dim_0    (mu_Y_dim_0) int64 4kB 0 1 2 3 4 5 ... 494 495 496 497 498 499
│         * mu_M_dim_0    (mu_M_dim_0) int64 4kB 0 1 2 3 4 5 ... 494 495 496 497 498 499
│       Data variables:
│           beta_Y        (chain, draw, Y_predictors) float64 96kB 0.01113 ... 0.3383
│           beta_M        (chain, draw, M_predictors) float64 64kB -0.04368 ... 0.4884
│           sigma_M       (chain, draw) float64 32kB 0.3083 0.2978 ... 0.2982 0.3008
│           sigma_Y       (chain, draw) float64 32kB 0.3004 0.311 0.3 ... 0.3008 0.303
│           mu_Y          (chain, draw, mu_Y_dim_0) float64 16MB 0.5574 ... -0.8974
│           mu_M          (chain, draw, mu_M_dim_0) float64 16MB 0.108 ... -0.7418
│       Attributes:
│           created_at:                 2026-09-30T20:04:11.935110+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.2.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.2.0
│           sample_dims:                ['chain', 'draw']
│           sampling_time:              2.6495630741119385
│           tuning_steps:               1000
├── Group: /sample_stats
│       Dimensions:                (chain: 4, draw: 1000)
│       Coordinates:
│         * chain                  (chain) int64 32B 0 1 2 3
│         * draw                   (draw) int64 8kB 0 1 2 3 4 5 ... 995 996 997 998 999
│       Data variables: (12/18)
│           max_energy_error       (chain, draw) float64 32kB -0.41 0.5627 ... -0.1698
│           tree_depth             (chain, draw) int64 32kB 2 3 3 3 3 3 ... 3 3 3 3 3 3
│           energy_error           (chain, draw) float64 32kB -0.1059 ... -0.168
│           diverging              (chain, draw) bool 4kB False False ... False False
│           perf_counter_start     (chain, draw) float64 32kB 1.083e+06 ... 1.083e+06
│           energy                 (chain, draw) float64 32kB 248.6 249.4 ... 246.7
│           ...                     ...
│           acceptance_rate        (chain, draw) float64 32kB 1.0 0.7703 ... 1.0 0.9995
│           divergences            (chain, draw) int64 32kB 0 0 0 0 0 0 ... 0 0 0 0 0 0
│           reached_max_treedepth  (chain, draw) bool 4kB False False ... False False
│           largest_eigval         (chain, draw) float64 32kB nan nan nan ... nan nan
│           index_in_trajectory    (chain, draw) int64 32kB 2 2 2 -7 -4 7 ... 4 1 -3 3 3
│           perf_counter_diff      (chain, draw) float64 32kB 7.879e-05 ... 0.0001337
│       Attributes:
│           created_at:                 2026-09-30T20:04:11.942767+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.2.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.2.0
│           sample_dims:                ['chain', 'draw']
│           sampling_time:              2.6495630741119385
│           tuning_steps:               1000
├── Group: /observed_data
│       Dimensions:  (M_dim_0: 500, Y_dim_0: 500)
│       Coordinates:
│         * M_dim_0  (M_dim_0) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
│         * Y_dim_0  (Y_dim_0) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
│       Data variables:
│           M        (M_dim_0) float64 4kB 0.5615 -0.2514 0.1594 ... -0.9562 -0.5004
│           Y        (Y_dim_0) float64 4kB 0.5228 -0.7319 0.2283 ... -1.149 -1.203
│       Attributes:
│           created_at:                 2026-09-30T20:04:11.945196+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.2.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.2.0
│           sample_dims:                []
├── Group: /constant_data
│       Dimensions:  (X_dim_0: 500)
│       Coordinates:
│         * X_dim_0  (X_dim_0) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
│       Data variables:
│           X        (X_dim_0) float64 4kB 0.3047 -1.04 0.7505 ... -0.3356 -1.991 -1.495
│       Attributes:
│           created_at:                 2026-09-30T20:04:11.945885+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.2.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.2.0
│           sample_dims:                []
└── Group: /log_likelihood
        Dimensions:  (chain: 4, draw: 1000, M_dim_0: 500, Y_dim_0: 500)
        Coordinates:
          * chain    (chain) int64 32B 0 1 2 3
          * draw     (draw) int64 8kB 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999
          * M_dim_0  (M_dim_0) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
          * Y_dim_0  (Y_dim_0) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
        Data variables:
            M        (chain, draw, M_dim_0) float64 16MB -0.8244 -0.247 ... -0.03983
            Y        (chain, draw, Y_dim_0) float64 16MB 0.2772 0.06512 ... -0.2345
        Attributes:
            created_at:                 2026-09-30T20:04:12.372087+00:00
            creation_library:           ArviZ
            creation_library_version:   1.2.0
            creation_library_language:  Python
            inference_library:          pymc
            inference_library_version:  6.2.0
            sample_dims:                ['chain', 'draw']
model.effects_summary()
mean sd hdi_3% hdi_97%
name
a 0.497194 0.014603 0.468207 0.523708
b 0.762801 0.046157 0.676966 0.851365
c 0.334834 0.026907 0.284274 0.385843

The posterior means should be close to the true values (a = 0.5, b = 0.8, c = 0.3).

Specifying priors up front

If you already know your priors before fitting, you can pass them directly to model() instead of using the two-step inspect-then-refine approach.

model_upfront = pathmc.model(
    spec,
    data=df,
    priors={
        "beta_M": Prior("Normal", mu=0, sigma=2),
        "beta_Y": Prior("Normal", mu=0, sigma=2),
        "sigma_M": Prior("HalfNormal", sigma=0.5),
        "sigma_Y": Prior("HalfNormal", sigma=0.5),
    },
)

model_upfront.equations()

\begin{aligned} \beta_{M} &\sim \text{Normal}(mu=0,\, sigma=2) \\ \sigma_{M} &\sim \text{HalfNormal}(sigma=0.5) \\ \beta_{Y} &\sim \text{Normal}(mu=0,\, sigma=2) \\ \sigma_{Y} &\sim \text{HalfNormal}(sigma=0.5) \\[6pt] \mu_{M} &= \beta_{0,\,M} + a \cdot \mathrm{X} \\ \mathrm{M} &\sim \text{Normal}(\mu_{M},\, \sigma_{M}) \\ \mu_{Y} &= \beta_{0,\,Y} + b \cdot \mathrm{M} + c \cdot \mathrm{X} \\ \mathrm{Y} &\sim \text{Normal}(\mu_{Y},\, \sigma_{Y}) \end{aligned}

Both approaches produce the same compiled model.

Using different distribution families

The Prior class supports any PyMC distribution. When domain knowledge suggests a specific shape — for example, that a coefficient should be positive, or bounded — you can express that directly.

model_informative = pathmc.model(
    spec,
    data=df,
    priors={
        "sigma_M": Prior("Exponential", lam=3),
        "sigma_Y": Prior("Exponential", lam=3),
    },
)

model_informative.equations()

\begin{aligned} \beta_{M} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{M} &\sim \text{Exponential}(lam=3) \\ \beta_{Y} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{Y} &\sim \text{Exponential}(lam=3) \\[6pt] \mu_{M} &= \beta_{0,\,M} + a \cdot \mathrm{X} \\ \mathrm{M} &\sim \text{Normal}(\mu_{M},\, \sigma_{M}) \\ \mu_{Y} &= \beta_{0,\,Y} + b \cdot \mathrm{M} + c \cdot \mathrm{X} \\ \mathrm{Y} &\sim \text{Normal}(\mu_{Y},\, \sigma_{Y}) \end{aligned}

TipDiscovering valid prior keys

Call .equations() on any model to see both the structural equations and the full list of prior parameter names you can customize. Use .equations(show="priors") to see only the priors. The keys depend on your model specification — a model with random slopes will have mu_slope_* and sigma_slope_* keys, while a model with transforms will include the transform parameter names.

Summary

  • Default priors (Normal(0, 10) for coefficients, HalfNormal(1) for scales) work out of the box and require zero configuration.
  • .equations() shows the structural equations and the current prior specification for every model parameter. Use show="priors" to see only the priors.
  • .sample_prior_predictive() generates data from the prior to verify that your assumptions produce plausible outcomes.
  • .set_priors() merges overrides into the current prior configuration and recompiles the model. Only the specified parameters change; everything else keeps its current value.
  • model(priors=...) accepts priors at creation time for a one-step workflow.
  • Any PyMC distribution can serve as a prior via Prior("DistName", ...) from pymc_extras.
  • The iterative workflow — inspect, check, refine, re-check, sample — is the principled approach to prior specification in Bayesian modeling. See Bayesian Workflow for the full end-to-end cycle from model specification through causal queries.
NoteAdapt this to your domain

Consider the quantities in your own causal model. What range of coefficient values would be scientifically absurd? If a 1-unit change in an input could never plausibly produce a 50-unit change in the outcome, a Normal(0, 10) prior is too generous.

  • Clinical trials: effect sizes above 1–2 standard deviations are rare; a Normal(0, 1) or Normal(0, 2) on standardized coefficients is often appropriate.
  • Marketing: elasticities (percent change in sales per percent change in spend) are typically in [0, 0.5]; a HalfNormal(sigma=0.3) keeps the prior realistic.
  • Economics: use prior studies or meta-analyses to set informative priors on known structural parameters.