Model-Based Data Simulation

Generate synthetic datasets from a pathmc model with known parameter values using pathmc.simulate().

Author

Benjamin Vincent

Most foundational examples in these docs build their datasets with explicit NumPy array code that anyone can audit. But pathmc can also generate data directly from a specified model, guaranteeing that the simulated dataset is perfectly consistent with the model’s generative structure.

This is useful for simulate-and-recover workflows: choose true parameter values, generate a dataset, fit the model, and verify that the posterior concentrates around truth. It is also a quick way to prototype a model before collecting real data.

The pathmc.simulate() function builds a generative PyMC model from the specification, fixes all parameters at user-provided values via pm.do(), and draws one simulated dataset.

MMM examples run on simulate()

The applied-models notebooks are moving to pathmc.simulate() for their ground-truth DGPs: MMM with Latent Brand Awareness simulates a latent AR(1) mediator via the latent= and families= arguments, and pairs with Media Mix Models, which walks through panel adstock/saturation structures under partial pooling. Everything shown on this page extends to those settings through the same params interface.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pathmc

FIG_WIDTH = 7
FIG_HEIGHT = 4

Simple regression

The simplest case: one exogenous predictor X and one outcome Y.

Figure 1: Single-equation regression: X causes Y.

Simulate

Define true parameter values and create the exogenous data. The parameter names match the PyMC model internals — beta_Y is the coefficient vector (ordered [Intercept, X]) and sigma_Y is the residual standard deviation.

seed = sum(map(ord, "simple regression"))
rng = np.random.default_rng(seed)

truth_simple = {
    "beta_Y": [2.0, 0.8],  # Intercept=2.0, slope=0.8
    "sigma_Y": 1.0,
}

n = 300
exog = pd.DataFrame({"X": rng.normal(size=n)})

df_simple = pathmc.simulate(
    "Y ~ X",
    data=exog,
    params=truth_simple,
    random_seed=rng,
)
df_simple.head()

	X	Y
0	-0.077591	2.910311
1	-0.704216	1.926878
2	-0.688257	3.130698
3	0.683797	1.672497
4	-1.683856	0.715676

The returned DataFrame contains the original exogenous column X plus the simulated outcome Y.

Code

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))
ax.scatter(df_simple["X"], df_simple["Y"], alpha=0.3, s=15, color="C0")
x_grid = np.linspace(df_simple["X"].min(), df_simple["X"].max(), 100)
y_true = truth_simple["beta_Y"][0] + truth_simple["beta_Y"][1] * x_grid
ax.plot(x_grid, y_true, "k--", lw=2, label="True: Y = 2.0 + 0.8X")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.legend()
plt.show()

Figure 2: Simulated data from Y ~ X with known intercept (2.0) and slope (0.8). Black dashed line shows the true regression function.

Parameter recovery

Fit the model on the simulated data and check that the posterior covers the true values.

model_simple = pathmc.model("Y ~ a*X", data=df_simple)
model_simple.fit(random_seed=42)
model_simple.effects_summary()

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta_Y, sigma_Y]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.

	mean	sd	hdi_3%	hdi_97%
name
a	0.8276	0.059604	0.716054	0.939373

The a coefficient should be close to 0.8. The full summary includes the intercept and residual scale:

model_simple.summary()

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
beta_Y[Intercept]	1.923114	0.061244	1.825989	2.019017	6554.115772	3358.488004	1.001088	0.000754	0.000552
beta_Y[X]	0.827600	0.059604	0.732422	0.921756	6015.486616	3221.375608	1.000979	0.000769	0.000566
sigma_Y	1.050665	0.042242	0.985317	1.120334	5879.244756	3013.500747	0.999858	0.000548	0.000401
mu_Y[0]	1.858900	0.061277	1.761407	1.955420	6567.258581	3386.499814	1.000817	0.000755	0.000556
mu_Y[1]	1.340305	0.073182	1.225129	1.455969	6428.835309	3327.645796	1.001188	0.000918	0.000682
...	...	...	...	...	...	...	...	...	...
mu_Y[295]	-0.313806	0.170596	-0.579586	-0.043726	6309.072293	3216.276711	1.000855	0.002156	0.001595
mu_Y[296]	1.105409	0.083659	0.974959	1.238437	6410.583403	3463.255555	1.001084	0.001051	0.000780
mu_Y[297]	2.045700	0.062143	1.945574	2.142614	6504.073572	3404.598401	1.000770	0.000767	0.000557
mu_Y[298]	2.329204	0.068670	2.215815	2.437176	6385.640422	3208.196625	1.001115	0.000858	0.000612
mu_Y[299]	2.367600	0.069966	2.251583	2.477707	6441.562083	3300.903849	1.001313	0.000876	0.000624

303 rows × 9 columns

Mediation chain

A two-equation model where X affects Y both directly and indirectly through a mediator M. Simulation must generate M first, then wire it into Y’s structural equation — pathmc.simulate() handles this automatically via the generative PyMC graph.

Figure 3: Mediation DAG: X affects Y directly and indirectly through M.

Simulate

truth_med = {
    "beta_M": [0.0, 0.5],  # M = 0.0 + 0.5*X + noise
    "sigma_M": 0.5,
    "beta_Y": [1.0, 0.8, 0.3],  # Y = 1.0 + 0.8*M + 0.3*X + noise
    "sigma_Y": 0.5,
}

# Derived truths
true_indirect = truth_med["beta_M"][1] * truth_med["beta_Y"][1]  # a*b
true_direct = truth_med["beta_Y"][2]  # c
true_total = true_direct + true_indirect  # c + a*b

print(f"True indirect effect (a×b): {true_indirect:.2f}")
print(f"True direct effect (c):     {true_direct:.2f}")
print(f"True total effect:          {true_total:.2f}")

True indirect effect (a×b): 0.40
True direct effect (c):     0.30
True total effect:          0.70

exog_med = pd.DataFrame({"X": rng.normal(size=500)})

df_med = pathmc.simulate(
    """
    M ~ X
    Y ~ M + X
    """,
    data=exog_med,
    params=truth_med,
    random_seed=rng,
)
df_med.head()

	X	M	Y
0	0.674129	0.890260	1.952226
1	0.433372	-0.091577	1.564749
2	0.773107	0.273978	1.851593
3	-0.469540	-0.431065	1.050672
4	-0.317956	-0.748322	0.671664

The correlation structure should reflect the causal chain — X and M are correlated (through a), and M and Y are correlated (through b and the X → M → Y path):

df_med[["X", "M", "Y"]].corr().round(2)

	X	M	Y
X	1.00	0.72	0.73
M	0.72	1.00	0.82
Y	0.73	0.82	1.00

Parameter recovery

spec_med = """
M ~ a*X
Y ~ b*M + c*X
indirect := a*b
total := c + a*b
"""

model_med = pathmc.model(spec_med, data=df_med)
idata_med = model_med.fit(random_seed=42)
model_med.effects_summary()

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 2 seconds.

	mean	sd	hdi_3%	hdi_97%
name
a	0.489955	0.021226	0.450047	0.529284
b	0.850396	0.047104	0.761491	0.936701
c	0.272720	0.032858	0.212078	0.333703
indirect	0.416674	0.029567	0.359805	0.471200
total	0.689395	0.028851	0.631791	0.740585

The indirect row should be close to {true_indirect} and total close to {true_total}. The path-specific effect API provides another way to check:

print(model_med.effect("X -> M -> Y"))
print(model_med.effect("X -> Y"))

EffectResult('X -> M -> Y', mean=0.4167, 94% HDI=[0.3598, 0.4712])
EffectResult('X -> Y', mean=0.2727, 94% HDI=[0.2121, 0.3337])

Code

import arviz as az

idata = idata_med
params_to_check = {
    "a (X→M)": ("beta_M", 1, truth_med["beta_M"][1]),
    "b (M→Y)": ("beta_Y", 1, truth_med["beta_Y"][1]),
    "c (X→Y)": ("beta_Y", 2, truth_med["beta_Y"][2]),
}

fig, axes = plt.subplots(1, 3, figsize=(FIG_WIDTH, FIG_HEIGHT * 0.8))
for ax, (label, (param, idx, true_val)) in zip(axes, params_to_check.items()):
    draws = (
        idata
        .posterior[param]
        .sel({
            f"{param.split('_')[1]}_predictors": idata
            .posterior[param]
            .coords[f"{param.split('_')[1]}_predictors"]
            .values[idx]
        })
        .values.flatten()
    )
    x_kde, y_kde, _ = az.kde(draws)
    ax.plot(x_kde, y_kde, color=None, lw=2)
    ax.fill_between(x_kde, y_kde, alpha=0.3, color=None)
    ax.axvline(true_val, color="k", ls="--", lw=1.5)
    ax.set_title(label)
    ax.set_yticks([])
fig.tight_layout()
plt.show()

Figure 4: Parameter recovery for the mediation model. Black dashed lines mark the true values used in simulation.

Binary outcomes

Simulation works with non-Gaussian families. For a Bernoulli outcome, the coefficients are on the logit scale — the same parameterization used internally by the generative model. No separate link-function bookkeeping is needed.

Figure 5: Binary outcome model: X1 and X2 jointly predict a binary response Y.

Simulate

With Bernoulli family, there is no sigma_Y — the only parameters are the regression coefficients. The intercept controls the base log-odds, and the slopes control how each predictor shifts the probability.

truth_bin = {
    "beta_Y": [-0.5, 1.2, -0.8],  # logit(P(Y=1)) = -0.5 + 1.2*X1 - 0.8*X2
}

exog_bin = pd.DataFrame({
    "X1": rng.normal(size=500),
    "X2": rng.normal(size=500),
})

df_bin = pathmc.simulate(
    "Y ~ X1 + X2",
    data=exog_bin,
    params=truth_bin,
    families={"Y": "bernoulli"},
    random_seed=rng,
)
df_bin.head()

	X1	X2
0	-1.275853	0.316749
1	0.663898	-0.353732
2	-0.725828	0.926820
3	-0.002459	-0.722710
4	-1.101098	0.975424

print(f"Outcome prevalence: {df_bin['Y'].mean():.2%}")
print(f"Outcome values: {sorted(df_bin['Y'].unique())}")

Outcome prevalence: 43.00%
Outcome values: [np.int64(0), np.int64(1)]

Code

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))
colors = ["C3" if y == 1 else "C0" for y in df_bin["Y"]]
ax.scatter(df_bin["X1"], df_bin["X2"], c=colors, alpha=0.3, s=15)

# Decision boundary: -0.5 + 1.2*X1 - 0.8*X2 = 0  →  X2 = (-0.5 + 1.2*X1) / 0.8
x1_grid = np.linspace(df_bin["X1"].min(), df_bin["X1"].max(), 100)
x2_boundary = (truth_bin["beta_Y"][0] + truth_bin["beta_Y"][1] * x1_grid) / -truth_bin[
    "beta_Y"
][2]
ax.plot(x1_grid, x2_boundary, "k--", lw=2, label="Decision boundary (P=0.5)")
ax.set_xlabel("X1")
ax.set_ylabel("X2")

from matplotlib.lines import Line2D

legend_elements = [
    Line2D(
        [0],
        [0],
        marker="o",
        color="w",
        markerfacecolor="C0",
        markersize=8,
        label="Y = 0",
    ),
    Line2D(
        [0],
        [0],
        marker="o",
        color="w",
        markerfacecolor="C3",
        markersize=8,
        label="Y = 1",
    ),
    Line2D([0], [0], color="k", ls="--", lw=2, label="Decision boundary"),
]
ax.legend(handles=legend_elements)
plt.show()

Figure 6: Simulated binary outcome colored by Y. The decision boundary (P(Y=1) = 0.5) is shown as a black dashed line.

Parameter recovery

model_bin = pathmc.model("Y ~ X1 + X2", data=df_bin, families={"Y": "bernoulli"})
model_bin.fit(random_seed=42)
model_bin.summary()

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta_Y]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
beta_Y[Intercept]	-0.422015	0.115591	-0.612535	-0.242165	4722.482116	3291.780316	1.000253	0.001682	0.001220
beta_Y[X1]	1.356547	0.146042	1.126517	1.593402	3746.799199	3204.280243	1.001526	0.002387	0.001754
beta_Y[X2]	-0.976825	0.127714	-1.189225	-0.779720	3951.930310	3044.602729	0.999942	0.002043	0.001419
mu_Y[0]	-2.462177	0.257041	-2.877148	-2.059632	3318.118857	3102.248808	1.002870	0.004463	0.003275
mu_Y[1]	0.824128	0.147049	0.590466	1.058438	5537.488159	3390.017850	1.000045	0.001984	0.001420
...	...	...	...	...	...	...	...	...	...
mu_Y[495]	-3.052298	0.307241	-3.540098	-2.571213	3241.180303	2939.245889	1.002767	0.005409	0.003963
mu_Y[496]	0.795477	0.153266	0.550873	1.040264	5697.441520	3384.646914	1.000499	0.002034	0.001439
mu_Y[497]	-0.770650	0.156159	-1.025557	-0.526624	3896.815931	2923.476153	0.999474	0.002501	0.001774
mu_Y[498]	-1.827077	0.199311	-2.145993	-1.516371	3264.541041	3017.562586	1.001807	0.003498	0.002572
mu_Y[499]	-0.109257	0.190052	-0.417387	0.190816	4667.800964	3239.552958	1.000442	0.002781	0.001916

503 rows × 9 columns

The beta_Y coefficients should recover the true logit-scale values: intercept ≈ −0.5, X1 ≈ 1.2, X2 ≈ −0.8.

Discovering parameter names

The parameter names passed to simulate() must match the PyMC variable names in the compiled model. To discover them for a new spec, build a temporary model with placeholder data and inspect its structure:

spec = """
M ~ X
Y ~ M + X
"""
placeholder = pd.DataFrame({"X": [0.0], "M": [0.0], "Y": [0.0]})
tmp = pathmc.model(spec, data=placeholder)
tmp.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} + \mathrm{X} \\ \mathrm{M} &\sim \text{Normal}(\mu_{M},\, \sigma_{M}) \\ \mu_{Y} &= \beta_{0,\,Y} + \mathrm{M} + \mathrm{X} \\ \mathrm{Y} &\sim \text{Normal}(\mu_{Y},\, \sigma_{Y}) \end{aligned}

`simulate_params_template()`

pathmc.simulate_params_template() replaces this manual step: it compiles the same zero-filled placeholder generative model that simulate() builds internally and returns every required parameter name with its expected shape and dtype — without drawing any data.

template = pathmc.simulate_params_template(
    """
    M ~ X
    Y ~ M + X
    """,
    data=pd.DataFrame({"X": [0.0]}),
)
for name, desc in sorted(template.items()):
    print(f"{name:>8}: {desc['kind']:>6} {str(desc['shape']):>5} ({desc['dtype']})")

  beta_M: vector  (2,) (float64)
  beta_Y: vector  (3,) (float64)
 sigma_M: scalar    () (float64)
 sigma_Y: scalar    () (float64)

Each value describes how to fill params: "kind" is "scalar", "vector", or "matrix", "shape" lists the axis lengths, and "dtype" gives the NumPy dtype. Panel models accept panel= and pooling= exactly as model() does; under pooling="partial", unit-indexed hierarchical parameters are reported with concrete shapes derived from the placeholder data (e.g. alpha_{var} has one entry per panel unit), and transform terms contribute their named hyperparameters (e.g. the adstock decay=theta_tv argument appears as theta_tv).

Summary

pathmc.simulate() generates data from a model specification with known parameter values, using the same generative PyMC graph that powers estimation and do() queries.
Parameters are fixed via pm.do() and data is drawn via pm.draw(), guaranteeing consistency with the model’s internal structure — including multi-equation chains and link functions.
The returned DataFrame contains the original exogenous columns plus simulated endogenous columns.
Supported for all families (gaussian, bernoulli, poisson, negbinomial, studentt) and multi-equation models.
Use equations() on a placeholder model to discover both the structural equations and expected parameter names.
pathmc.simulate_params_template() lists every parameter name, shape, and dtype a spec requires — including hierarchical, transform, residual-covariance, and hsgp parameters — so params can be authored without trial and error.

When to use numpy simulation instead

pathmc.simulate() is convenient for simulate-and-recover and prototyping. However, for pedagogical datasets where readers need to see and audit the exact data-generating process, explicit NumPy code is often clearer — every line maps directly to a structural equation, and there is no abstraction layer to look through. The two approaches should produce equivalent data (up to random variation) when the DGP matches the model spec.