Simulating Data with PyMC

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

This walkthrough draws on Ricardo Vieira and Tomás Capretto’s PyMC Labs simulation tutorial. The snippets below target PyMC 5.24.0 APIs and require NumPy, pandas, SciPy, seaborn, and Matplotlib. The diamonds example uses seaborn’s public example data; loading it requires network access or its cached copy.

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:

import numpy as np
import pandas as pd
import pymc as pm
import scipy.stats as st
import seaborn as sns
import matplotlib.pyplot as plt
from pymc.pytensorf import compile as compile_random
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.

Are PyMC distributions 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 returns a three-way split summing to one. Given the specified weights, the first row usually places most weight on its third component and the second on its first. Exact draws depend on the software stack.

Meta-distributions: truncation and mixtures

PyMC’s Truncated API supports specified univariate distributions subject to its API constraints:

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);

How do mixture distributions combine 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

Each selected draw shares the sampled categorical index and therefore uses the corresponding Normal mean. This graph dependence matters more than a particular seeded result.

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)

The event-size sum changes with the sampled Poisson count. Seeded values are not portable output guarantees.

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.

df = sns.load_dataset("diamonds").dropna(subset=["cut", "price"])
cut_labels = pd.Categorical(df["cut"])
cut_idxs = cut_labels.codes
n_cuts = len(cut_labels.categories)
model_coords = {
    "obs": range(len(df)),
    "components": (0, 1, 2),
    "cut": list(cut_labels.categories),
}
with pm.Model(coords=model_coords) as diamond_model:
    prior_dims = ("cut", "components")
    weight_prior = pm.Dirichlet("mix_weights", np.ones((n_cuts, 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 fitted dictionary has one row per category and one column per component. Inspect labels rather than assuming component identities match across cuts or runs.

These parameter summaries depend on the fit and category ordering; no exact numerical output is guaranteed.

Validating the fit against real data

PyMC’s do transformation returns a distinct model. Retrieve its cloned outcome rather than drawing the original graph variable:

fitted_model = pm.do(diamond_model, fit)
draws = pm.draw(fitted_model["price"], random_seed=10, draws=20)

Here each draw has one simulated price per observation, with the observation’s cut index. Compare draws[:, cut_idxs == i].ravel() with observed prices in category i. Histogram and ECDF agreement checks fitted marginals; it is not a holdout test, joint-distribution validation, or causal validation. Check a new reference dataset and revise the model when those checks fail.

Performance notes for draw-heavy workflows

For repeated draws from the same graph, compile the random function once and reuse it. This avoids repeatedly compiling through separate draw calls. Measure performance for the actual shape and backend.

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

The two calls advance the compiled random stream; exact values can differ across library versions.

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_random(inputs=[], outputs=draw_dist, random_seed=11)
draw_fn()

The call returns an array with the requested two-element shape.

Why does the fit-then-check step matter 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 complete primary tutorial, versioned do API, and compilation API provide implementation detail. Discuss validation of a decision model when marginal agreement does not answer the intended business question.