Skip to main content

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 with nu_{var}. Scan lag-carry reconstructs each state as mu + sigma_{var} * e, where e is a unit normal innovation, so the innovation scale is sigma_{var} and nu is not used. For residual covariances (~~), supply "chol_{block_name}" as the packed lower-triangular Cholesky vector of length k(k+1)/2 (PyMC LKJCholeskyCov packing); 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 raise NotImplementedError. Hierarchical unit effects (alpha_{var}, slope_{lhs}_{var}) are passed as the centered unit vectors, not the *_offset free parameters the sampler uses. For hsgp() terms, supply "ell_{lhs}_{var}", "eta_{lhs}_{var}", and "beta_hsgp_{lhs}_{var}" (length m). 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. Use pathmc.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 sigma parameter; 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 uses lag() 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 member sigma parameter 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 the alpha_{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-cell beta_{var} vectors and the dim-indexed mu_{var}_{dim} / scalar sigma_{var}_{dim} hyperpriors), or None for 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 but panel= 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’]