Custom Priors
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.
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()
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()
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}
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. Useshow="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", ...)frompymc_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.
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)orNormal(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.