Simulating Data with PyMC

Edited by Avi Yashchin, CEO and Co-Founder · August 20, 2026

A simulated output is only as trustworthy as the process that generated it. Before a data science or analytics leader commits budget to a decision built on a simulated experiment, the underlying data-generating process needs to be specified explicitly, fit to real data, and checked against that data. PyMC is a probabilistic programming library built for exactly that workflow, and its simulation tools show the discipline in miniature.

Draws from a specified distribution

PyMC models a highly structured data-generating process, useful even outside formal Bayesian inference: for simulation in optimization routines, risk analysis, and research design.

A distribution is defined with explicit parameters, then sampled with draw:

spend = pm.Gamma.dist(alpha=2, beta=1)
spend_draws = pm.draw(spend, draws=1000, random_seed=1)
sns.histplot(spend_draws);

The same distribution, reparametrized

The same Gamma distribution can be defined with a mean/standard-deviation parametrization instead of shape/rate:

spend_alt = pm.Gamma.dist(mu=2, sigma=2**0.5)
spend_alt_draws = pm.draw(spend_alt, draws=1000, random_seed=2)
sns.histplot(spend_alt_draws);

PyMC converts between equivalent parametrizations automatically, so the choice is about what's easiest to reason about, not a modeling compromise.

Vectorized by default

Not every scientific-computing distribution allows array-style broadcasting across its parameters. Every distribution in the PyMC distributions API does. A two-row Dirichlet call in SciPy raises an error:

cut_weights = [[1, 5, 100], [100, 5, 1]]
try:
    st.dirichlet(cut_weights).rvs()
except ValueError as broadcast_error:
    print(broadcast_error)

SciPy raises: parameter vector a must be one dimensional, but its shape is (2, 3).

Passing the array directly

PyMC distributions are vectorized, so the same call runs directly:

rows = pm.Dirichlet.dist(cut_weights)
pm.draw(rows, random_seed=3)

Each row draws its own three-way split that sums to 1: the first row (weights skewed toward 100) lands near [0.0017, 0.0396, 0.9586], and the second (weights mirrored) lands near [0.9500, 0.0418, 0.0082].

Meta-distributions: truncation and mixtures

Real-world quantities are often bounded, or drawn from more than one underlying process. PyMC's Truncated class can bound any univariate distribution:

unbounded = pm.Lognormal.dist(0, 1)
bounded = pm.Truncated.dist(unbounded, upper=3)
bounded_draws = pm.draw(bounded, draws=10_000, random_seed=4)
sns.histplot(bounded_draws);

Mixing multiple components

Mixture distributions combine multiple components with weights. Combining two Normals with weights [0.3, 0.7] draws roughly 30% of samples from the first component and 70% from the second:

component_low = pm.Normal.dist(-1, 1)
component_high = pm.Normal.dist(1, 0.5)
mix_weights = [0.3, 0.7]
mix = pm.Mixture.dist(w=mix_weights, comp_dists=[component_low, component_high])
mix_draws = pm.draw(mix, draws=10_000, random_seed=5)
sns.histplot(mix_draws);

Mixtures as a starting state

The same mixture machinery composes with random walks, letting a mixture serve as the initial state for a multi-step process:

start_weights = [0.3, 0.7]
start_low = pm.Beta.dist(1, 1)
start_high = pm.Normal.dist(100, 0.5)
start_dist = pm.Mixture.dist(w=start_weights, comp_dists=[start_low, start_high])

step_dist = pm.StudentT.dist(mu=0, sigma=1, nu=4)
walk = pm.RandomWalk.dist(init_dist=start_dist, innovation_dist=step_dist, steps=1000)

walk_draws = pm.draw(walk, draws=5, random_seed=6)
for path_draw in walk_draws:
    plt.plot(path_draw)
plt.xlabel("t");

Dependent variables in one draw

Simulated variables are often not independent of each other. This example samples a categorical index, then uses it to select from a vector of three Normals with means [-100, 0, 100]:

category_idx = pm.Categorical.dist(p=[.1, .3, .6])
selected = pm.Normal.dist(mu=[-100, 0, 100], sigma=1)[category_idx]
idx_samples, selected_samples = pm.draw([category_idx, selected], draws=5, random_seed=7)
idx_samples, selected_samples

The five category draws come back as [2, 0, 1, 1, 0], and the matching selected values land near each drawn category's mean: around 99.6, -97.9, -1.9, -0.5, and -101.1.

One draw sizing another

A Poisson draw can also determine the shape of a downstream draw: how many Gamma-distributed events to sum.

event_count = pm.Poisson.dist(5)
event_sizes = pm.Gamma.dist(mu=10, sigma=2, shape=event_count)
pm.draw([event_count, event_sizes.sum()], draws=3, random_seed=8)

Across the three draws, the event counts come back [4, 10, 5], and the matching event-size sums come back roughly [36.4, 100.0, 57.8]: more Poisson events push the summed Gamma total higher.

When the process itself is unknown, infer it

The examples above start from a distribution the analyst already chose. A harder and more common case: a rough sense of what the data looks like, without a clear specification for how to simulate it. That is the setup behind fitting a model to real data, where the goal is realistic covariates whose marginals match an observed dataset.

Because PyMC also performs inference, the same framework can recover the parameters of an assumed structure from real data. Given price data across five diamond cuts (53,940 rows), a reasonable guess is a mixture of three LogNormal distributions per cut: 5 × 3 = 15 means, 15 standard deviations, and 15 mixture weights.

model_coords = {
    "obs": range(len(df)),
    "components": (0, 1, 2),
    "cut": cut_labels.codes,
}
with pm.Model(coords=model_coords) as diamond_model:
    prior_dims = ("cut", "components")
    weight_prior = pm.Dirichlet("mix_weights", np.ones((5, 3)), dims=prior_dims)
    cut_means = [7, 8, 9]
    mean_prior = pm.Normal("mix_means", dims=prior_dims, mu=cut_means, sigma=3)
    std_prior = pm.HalfNormal("mix_stds", sigma=2, dims=prior_dims)

    price = pm.Mixture(
        "price",
        w=weight_prior[cut_idxs],
        comp_dists=pm.LogNormal.dist(mu=mean_prior[cut_idxs], sigma=std_prior[cut_idxs]),
        observed=df["price"],
        dims="obs",
    )

A MAP point estimate

Fitting with find_MAP returns one local-optimum point estimate of weights, means, and standard deviations for each cut, with no uncertainty and no guarantee that component 1 means the same thing across cuts or runs:

with diamond_model:
    fit = pm.find_MAP(include_transformed=False)
fit

The dictionary returns one row of three means, three mixture weights, and three standard deviations per diamond cut. For the first cut alone, fit["mix_means"][0] is [6.66074039, 7.75610166, 9.05369536], fit["mix_weights"][0] is [0.35500674, 0.43518759, 0.20980567], and fit["mix_stds"][0] is [0.29090793, 0.58885309, 0.42603866]. The remaining four cuts return their own three-component fits in the same structure: five cuts times three components, each with its own recovered mean, weight, and standard deviation.

One local-optimum MAP point estimate for the first diamond cut, from fit["mix_means"][0], fit["mix_weights"][0], and fit["mix_stds"][0] above.

Validating the fit against real data

The next step is drawing from the model conditioned on the fitted MAP values and checking whether the simulated output resembles the real data it was fit to:

with pm.do(diamond_model, fit):
    draws = pm.draw(price, random_seed=10, draws=10_000)

Two checks make that comparison concrete: overlaying simulated and observed histograms per category, and overlaying their empirical CDFs. The marginal simulated histograms resembled the original dataset's, and the two ECDF lines tracked closely enough to trust the fit. When they don't line up, the answer is to revisit the model: more mixture components, or a different distribution family. Don't ship the simulated output anyway.

Performance notes for draw-heavy workflows

draw compiles a function, seeds it, and calls it in a Python loop. That's convenient for exploration, but costly in a loop, because it recompiles the same function on every call. For performance-critical code, compiling the underlying random function once with compile_pymc and reusing it avoids that cost:

draw_dist = pm.Normal.dist()
draw_fn = compile_pymc(inputs=[], outputs=draw_dist, random_seed=11)
draw_fn(), draw_fn()

(array(-0.10568235), array(-0.6541993))

Vectorizing the compiled call

Defining the distribution with its final shape and calling the compiled function once, instead of looping, lets the random number generation vectorize:

draw_dist = pm.Normal.dist(shape=(2,))
draw_fn = compile_pymc(inputs=[], outputs=draw_dist, random_seed=11)
draw_fn()

array([-0.10568235, -0.6541993 ])

Why the fit-then-check step matters for a decision

Everything above is general-purpose PyMC, not a benchmark or case study. But the discipline it demonstrates generalizes past any one tool: specify a structured process, fit it to real marginal data, then compare the simulated draws back against the observed distribution before treating the result as trustworthy.

Matching marginals like this establishes distributional realism for simulation inputs; it does not by itself validate the joint dependence structure or any causal effect estimated on the simulated population. The same standard causal experimentation applies to a simulated market's causal estimates, but those require their own separate validation before informing a launch, price, or message decision. A simulated experiment that was never checked against a real distribution is a model artifact, not evidence. Ship the wrong price or message on it, and the mistake surfaces only after the spend is committed.

The full PyMC notebook and API documentation cover the mechanics in more depth. The discipline that carries over is simpler than the code: fit the parameters, then validate against the observed data before the output drives a decision.