modeling.simulate()
Simulate data from a pathmc model with known parameter values.
Usage
modeling.simulate(
spec_string,
data,
params,
families=None,
latent=None,
panel=None,
pooling=None,
scaling=None,
random_seed=None
)Builds a generative PyMC model from the specification, fixes all parameter random variables at the values in params using pm.do(), and draws one simulated dataset.
This is useful for:
- Simulate-and-recover workflows: generate data from known parameters, build and fit the model, and verify that the posterior concentrates around the truth.
- Teaching: create pedagogical datasets with exact causal structure matching the model DAG.
- Power analysis: generate data under hypothesized effect sizes and check whether the model can detect them.
Parameters
spec_string: str-
Model specification in the pathmc DSL. All regressions define the generative structure (e.g.
"M ~ X\nY ~ M + X"). data: IntoFrame-
pandas or polars DataFrame containing the exogenous variables (predictors). Endogenous (outcome) columns need not be present — they will be simulated. If present, they are ignored.
params: dict[str, Any]-
True parameter values keyed by PyMC variable name. Typical keys are
"beta_{var}"(coefficient vector) and"sigma_{var}". For a Gaussian outcome,sigma_{var}is the residual standard deviation, which equals the distribution scale. For StudentT,sigma_{var}is the scale, not the residual standard deviation. simulate() leaves that supplied value unchanged so a posterior round-trip is not converted twice. Cross-sectional draws use that scale withnu_{var}. Scan lag-carry reconstructs each state asmu + sigma_{var} * e, whereeis a unit normal innovation, so the innovation scale issigma_{var}andnuis not used. For residual covariances (~~), supply"chol_{block_name}"as the packed lower-triangular Cholesky vector of lengthk(k+1)/2(PyMCLKJCholeskyCovpacking); block members’ simulated columns share the correlated residuals. Descendants of block members receive the realized noisy draws. Within-block directed edges (one member regressing on another) still raiseNotImplementedError. Hierarchical unit effects (alpha_{var},slope_{lhs}_{var}) are passed as the centered unit vectors, not the*_offsetfree parameters the sampler uses. Forhsgp()terms, supply"ell_{lhs}_{var}","eta_{lhs}_{var}", and"beta_hsgp_{lhs}_{var}"(lengthm). Transform parameters are keyed by their user-chosen DSL names, e.g.adstock(tv, decay=theta_tv)expects"theta_tv", and nested chains honor each link’s parameter. Usepathmc.model(...).equations()on a dummy dataset to discover expected names and shapes. families: dict[str, str] | None = None-
Per-variable distribution families (default
"gaussian"). Supports the same families as model():"gaussian","bernoulli","poisson","negbinomial","studentt","latent_normal". latent: list[str] | set[str] | None = None-
Variables to treat as latent (unobserved). These are compiled without a likelihood and their simulated values are included in the output. Deterministic latent nodes have no
sigmaparameter; stochastic latent nodes (families={"M": "latent_normal"}) do. panel: dict[str, str] | None = None-
Panel structure
{"unit": ..., "time": ...}mapping to column names, activating panel mode. Required when the spec useslag()terms (an error is raised otherwise, mirroring model()). With panel, temporal state starts cold: lagged endogenous terms and adstock carry begin at zero at each unit’s first time step, matching the compiler’s init semantics. Residual blocks (~~) are supported for complete observed Gaussian scan panels, with one globally shared covariance and independent innovations across units and time before propagation. Packed factors must be finite vectors with positive diagonals; no separate membersigmaparameter is required. Latent nodes and non-Gaussian equations cannot be combined with scan +~~. Supplied endogenous columns are ignored, including NaN columns. Convolution adstock compiled without panel treats row order as one series and warns. pooling: str | dict | None = None-
Pooling configuration forwarded to the compiler, identical to model():
"partial"for random intercepts per unit (then params must supply thealpha_{var}unit vectors, and their hierarchical means/scales even though they are clamped), a dict with"by_var"entries (then params must supply the per-cellbeta_{var}vectors and the dim-indexedmu_{var}_{dim}/ scalarsigma_{var}_{dim}hyperpriors), orNonefor complete pooling. scaling: pathmc.Scaling | ScalingFactors | None = None-
Inverse-scaling hook mirroring the
scaling=argument of model(). When given, exogenous channel columns are divided by the fitted factors before compilation — so params are interpreted in the same scaled units estimation uses — and every generated endogenous column is multiplied back by its factor, so outputs land in business units. Omit scaling entirely (the default) when your truth parameters already live in raw data units: nothing is scaled or inverted. Target specs must be grid-based ("fixed"/"divide") here because outcome values do not exist before simulation; for the exact estimation-time scales of a"max"/"mean"fit, pass the fitted model’s fitted_scaling object directly.
Returns
IntoFrame- Copy of data (same backend as the input) with simulated endogenous columns appended (including latent variables). Internal composite unit keys derived for multi-dimensional panels are not returned.
Raises
ValueError-
If required parameter values are missing from params, if a parameter array has the wrong shape, or if the spec contains
lag()terms butpanel=was omitted. NotImplementedError-
If a directed edge lies inside a
~~block. Drop the within-block path or the~~clause.
Examples
>>> import numpy as np
>>> import pandas as pd
>>> import pathmc
>>> rng = np.random.default_rng(42)
>>> exog = pd.DataFrame({"X": rng.normal(size=100)})
>>> df = pathmc.simulate(
... "Y ~ X",
... data=exog,
... params={"beta_Y": [2.0, 0.5], "sigma_Y": 1.0},
... random_seed=42,
... )
>>> list(df.columns)[‘X’, ‘Y’]