Monte Carlo Simulation

The Limits of Analytical Solutions

Many probability questions have clean analytical answers. The probability that a binomial random variable exceeds a threshold. The expected value of a sum of random variables. The confidence interval for a sample mean. These have formulas because the underlying distributions are well-understood and the questions are simple.

But real product metrics are rarely simple. Consider:

• What is the probability that total revenue next quarter exceeds $1M, given uncertainty in traffic, conversion rate, and average order value?
• What is the distribution of time to 1,000 signups, given a weekly growth rate that itself is uncertain?
• What is the probability that a new feature increases retention by at least 2 percentage points, given uncertainty in both the treatment and control groups?

These questions involve multiple uncertain quantities that interact nonlinearly. No closed-form solution exists. You could make simplifying assumptions (independence, normality, linearity) to force an analytical answer, but the assumptions would dominate the result and produce overconfident conclusions.

The Monte Carlo Alternative

Monte Carlo simulation replaces analytical derivation with computational sampling:

1. Model the uncertainty. Specify a distribution for each uncertain input (traffic ~ Normal(50000, 5000), conversion ~ Beta(50, 950), etc.).

2. Sample. Draw random values from each distribution.

3. Compute. Plug the sampled values into your metric formula.

4. Repeat. Do this thousands of times.

5. Analyze. The distribution of the computed metric across all simulations is your answer.

The output is not a single number — it is a full distribution. You can read off any quantile, any probability, any risk measure directly from the simulated samples.

A Simple Example: Revenue Uncertainty

Suppose monthly revenue depends on three uncertain inputs:

• Traffic: Normal(100,000, 10,000)
• Conversion rate: Beta(50, 950) (mean 5%)
• Average order value: Normal(80, 15)

Revenue = Traffic * Conversion * AOV. The distribution of this product has no simple closed form. Monte Carlo handles it directly:

import numpy as np

np.random.seed(42)
n_sims = 100_000

traffic = np.random.normal(100_000, 10_000, n_sims)
conversion = np.random.beta(50, 950, n_sims)
aov = np.random.normal(80, 15, n_sims)

revenue = traffic * conversion * aov

print(f'Mean revenue: ${revenue.mean():,.0f}')
print(f'5th percentile: ${np.percentile(revenue, 5):,.0f}')
print(f'95th percentile: ${np.percentile(revenue, 95):,.0f}')
print(f'P(revenue > $450K): {(revenue > 450_000).mean():.3f}')

The simulation gives you the full distribution: mean, spread, tail probabilities, and any other quantity you need. No formulas required.

When to Use Simulation

Use Monte Carlo when:

• The metric is a nonlinear function of uncertain inputs.
• The inputs have complex distributions (mixtures, truncated, empirical).
• You need tail probabilities that analytical approximations handle poorly.
• The problem is too complex for closed-form derivation.

Do not use Monte Carlo when:

• An analytical solution exists and is tractable. It is faster and exact.
• You need high precision in the far tails (beyond the 1st or 99th percentile). Simulation requires enormous sample sizes for rare events.
• The inputs are deterministic. Simulation adds noise without adding insight.

The art is knowing which situations fall into which category.

The Simulation Loop

Every Monte Carlo simulation follows the same structure:

import numpy as np

def run_simulation(n_sims, seed=42):
    rng = np.random.default_rng(seed)
    results = np.empty(n_sims)
    
    for i in range(n_sims):
        # 1. Sample inputs
        # 2. Compute metric
        # 3. Store result
        results[i] = compute_metric(rng)
    
    return results

The loop is conceptually simple, but the implementation matters. Vectorized operations (sampling all inputs at once) are orders of magnitude faster than Python loops:

# Vectorized: sample all at once
traffic = rng.normal(100_000, 10_000, n_sims)
conversion = rng.beta(50, 950, n_sims)
revenue = traffic * conversion * aov

# vs. loop: sample one at a time (100x slower)
for i in range(n_sims):
    t = rng.normal(100_000, 10_000)
    c = rng.beta(50, 950)
    revenue[i] = t * c * aov

Always vectorize when the metric computation allows it.

Choosing Input Distributions

The quality of a Monte Carlo simulation depends entirely on the input distributions. Garbage in, garbage out.

Strategies for specifying inputs:

1. Fit to historical data. Use the empirical distribution or fit a parametric distribution (normal, log-normal, gamma) to observed values.

2. Use domain knowledge. If traffic has grown 10% ± 3% per quarter historically, model next quarter's growth as Normal(0.10, 0.03).

3. Use Bayesian posteriors. If you have a posterior distribution from a previous analysis, sample from it directly.

4. Use empirical bootstrapping. Resample historical data with replacement when you do not want to assume a parametric form.

import numpy as np
from scipy import stats

# Fit a log-normal to historical revenue data
historical_revenue = np.array([...])  # your data
shape, loc, scale = stats.lognorm.fit(historical_revenue, floc=0)

# Sample from the fitted distribution
samples = stats.lognorm.rvs(shape, loc=loc, scale=scale, size=100_000)

Correlated Inputs

Real inputs are often correlated. Traffic and conversion may both spike during marketing campaigns. Revenue and costs move together. Ignoring correlation produces simulations that are too narrow (underestimating risk) or too wide (overestimating risk), depending on the sign of the correlation.

Use a multivariate distribution or a copula to model dependence:

import numpy as np

mean = [100_000, 0.05]
cov = [[10_000**2, 0.0002],
       [0.0002, 0.01**2]]

samples = np.random.multivariate_normal(mean, cov, size=100_000)
traffic, conversion = samples[:, 0], samples[:, 1]

The covariance matrix encodes both the individual variances and the correlation between inputs. The correlation coefficient ρ = Cov(X,Y) / (σ_X * σ_Y) should be estimated from historical data.

Convergence: How Many Simulations?

The precision of a Monte Carlo estimate improves with the square root of the number of simulations. The standard error of a probability estimate p̂ is:

SE = sqrt(p̂(1 - p̂) / n)

For p = 0.5 and n = 10,000, SE = 0.005. For n = 100,000, SE = 0.0016. To get one additional decimal digit of precision, you need 100x more simulations.

Practical guidance:

• 10,000 simulations for quick exploration and rough estimates.
• 100,000 simulations for production forecasts and decision-making.
• 1,000,000+ simulations for tail probabilities (beyond the 5th/95th percentile).

Check convergence by running the simulation with different sample sizes and comparing results:

import numpy as np

for n in [1_000, 10_000, 100_000]:
    samples = np.random.normal(100, 15, n)
    print(f'n={n}: mean={samples.mean():.2f}, std={samples.std():.2f}')

If the estimates stabilize as n increases, the simulation has converged. If they keep moving, you need more samples.

From Point Forecasts to Distributional Forecasts

Traditional forecasting produces a point estimate: "Next quarter's revenue will be $2.4M." This number is misleading because it implies certainty. The real question is: what is the distribution of possible outcomes?

Monte Carlo simulation produces a distributional forecast: "Next quarter's revenue has a 10% chance of being below $2.1M, a 50% chance of being below $2.4M, and a 10% chance of exceeding $2.7M."

This is more honest and more useful. It tells decision-makers the range of outcomes and the probability of each, enabling risk-aware planning.

A Quarterly Revenue Forecast

Suppose quarterly revenue depends on:

• New customers: Normal(5,000, 500)
• Revenue per customer: Log-normal(μ=4.5, σ=0.3) (median ~$90)
• Churn rate: Beta(20, 80) (mean 20%)

Revenue = New customers * Revenue per customer * (1 - Churn rate)

import numpy as np

np.random.seed(42)
n_sims = 100_000

new_customers = np.random.normal(5_000, 500, n_sims)
rev_per_customer = np.random.lognormal(4.5, 0.3, n_sims)
churn = np.random.beta(20, 80, n_sims)

revenue = new_customers * rev_per_customer * (1 - churn)

print(f'Median revenue: ${np.median(revenue):,.0f}')
print(f'10th percentile: ${np.percentile(revenue, 10):,.0f}')
print(f'90th percentile: ${np.percentile(revenue, 90):,.0f}')
print(f'P(revenue > $500K): {(revenue > 500_000).mean():.3f}')
print(f'P(revenue < $350K): {(revenue < 350_000).mean():.3f}')

The output is a full risk profile. You can answer any question about the forecast directly from the simulated distribution.

Scenario Analysis

Simulation enables scenario analysis: what happens to the forecast under different assumptions?

import numpy as np

def forecast_revenue(n_sims, mean_customers, std_customers, churn_alpha, churn_beta):
    new_customers = np.random.normal(mean_customers, std_customers, n_sims)
    rev_per_customer = np.random.lognormal(4.5, 0.3, n_sims)
    churn = np.random.beta(churn_alpha, churn_beta, n_sims)
    return new_customers * rev_per_customer * (1 - churn)

# Base case
base = forecast_revenue(100_000, 5_000, 500, 20, 80)

# Optimistic: more customers, lower churn
optimistic = forecast_revenue(100_000, 6_000, 500, 15, 85)

# Pessimistic: fewer customers, higher churn
pessimistic = forecast_revenue(100_000, 4_000, 500, 25, 75)

for name, dist in [('Base', base), ('Optimistic', optimistic), ('Pessimistic', pessimistic)]:
    print(f'{name}: median=${np.median(dist):,.0f}, '
          f'P(<$350K)={(dist < 350_000).mean():.3f}')

Scenario analysis shows how sensitive the forecast is to each assumption. If the pessimistic scenario has a high probability of missing the revenue target, you know which assumption to stress-test.

Time-Series Simulation

For metrics that evolve over time, simulate the entire trajectory:

import numpy as np

def simulate_dau_trajectory(n_days, initial_dau, growth_rate, growth_volatility, n_sims=10_000):
    rng = np.random.default_rng(42)
    trajectories = np.zeros((n_sims, n_days))
    trajectories[:, 0] = initial_dau
    
    for day in range(1, n_days):
        daily_growth = rng.normal(growth_rate, growth_volatility, n_sims)
        trajectories[:, day] = trajectories[:, day - 1] * (1 + daily_growth)
    
    return trajectories

trajectories = simulate_dau_trajectory(90, 50_000, 0.005, 0.02)

# Distribution of DAU at day 90
day_90_dau = trajectories[:, -1]
print(f'Median DAU at day 90: {np.median(day_90_dau):,.0f}')
print(f'5th-95th percentile: [{np.percentile(day_90_dau, 5):,.0f}, {np.percentile(day_90_dau, 95):,.0f}]')
print(f'P(DAU > 80K at day 90): {(day_90_dau > 80_000).mean():.3f}')

The trajectory simulation captures compounding uncertainty. Small daily variations accumulate into large uncertainty at the horizon. This is why long-term forecasts are inherently uncertain — not because the model is wrong, but because uncertainty compounds.

Problems Without Closed-Form Solutions

Some probability questions are easy to state but hard to solve analytically:

• What is the probability that variant B beats variant A by at least 2 percentage points in conversion rate?
• What is the probability that a user churns within 30 days, given a time-varying churn hazard?
• What is the probability that total infrastructure cost exceeds the budget, given uncertain traffic growth and per-request cost?

These questions involve multiple random variables, nonlinear relationships, and sometimes path dependence. Monte Carlo simulation handles them uniformly: model the inputs, compute the metric, count the fraction of simulations where the event occurs.

Example: Probability of Beating a Target

Suppose you want to know the probability that a new onboarding flow increases activation by at least 3 percentage points. You have posterior distributions from an A/B test:

• Control activation: Beta(40, 960) (mean 4%)
• Treatment activation: Beta(55, 945) (mean 5.5%)

import numpy as np

np.random.seed(42)
n_sims = 100_000

control = np.random.beta(40, 960, n_sims)
treatment = np.random.beta(55, 945, n_sims)

difference = treatment - control

print(f'P(treatment > control): {(difference > 0).mean():.3f}')
print(f'P(difference > 3pp): {(difference > 0.03).mean():.3f}')
print(f'Expected difference: {difference.mean():.4f}')
print(f'95% CI for difference: [{np.percentile(difference, 2.5):.4f}, {np.percentile(difference, 97.5):.4f}]')

The simulation answers questions that would require complex analytical work: the probability of any improvement, the probability of a clinically meaningful improvement, and the full distribution of the effect size.

Example: Time to Reach a Milestone

How many days until you reach 10,000 signups, given uncertain daily signup rates?

import numpy as np

def days_to_milestone(milestone, mean_daily, std_daily, n_sims=10_000):
    rng = np.random.default_rng(42)
    days = np.zeros(n_sims, dtype=int)
    
    for i in range(n_sims):
        total = 0
        day = 0
        while total < milestone:
            daily = max(0, rng.normal(mean_daily, std_daily))
            total += daily
            day += 1
        days[i] = day
    
    return days

days = days_to_milestone(10_000, 300, 50)

print(f'Median days to 10K: {np.median(days):.0f}')
print(f'90% CI: [{np.percentile(days, 5):.0f}, {np.percentile(days, 95):.0f}]')
print(f'P(within 40 days): {(days <= 40).mean():.3f}')

The while loop simulates each trajectory day by day until the milestone is reached. The distribution of completion times tells you the probability of hitting a launch deadline.

Variance Reduction Techniques

When you need precise estimates of rare events (e.g., P(metric > 5 sigma)), naive simulation requires millions of samples because the event occurs in only a few simulations. Variance reduction techniques improve efficiency:

1. Importance sampling. Sample from a distribution that overweights the rare event, then reweight the results. This concentrates samples in the region of interest.

2. Stratified sampling. Divide the input space into strata and sample proportionally from each. This ensures coverage of all regions.

3. Antithetic variates. For each random draw u, also use 1-u. This induces negative correlation between pairs, reducing variance.

import numpy as np

# Antithetic variates for a normal mean estimation
n = 50_000
u = np.random.random(n)
samples = np.concatenate([stats.norm.ppf(u), stats.norm.ppf(1 - u)])
print(f'Antithetic mean: {samples.mean():.6f}')

# vs. naive
naive = np.random.normal(0, 1, 2 * n)
print(f'Naive mean: {naive.mean():.6f}')

The antithetic estimator has lower variance because the positive and negative deviations cancel more effectively.

Validating Simulation Results

Always validate your simulation before trusting its output:

1. Check convergence. Run with increasing sample sizes and confirm the estimates stabilize.

2. Test with known answers. If you can solve a simplified version analytically, verify the simulation matches.

3. Inspect the input distributions. Plot histograms of the sampled inputs to confirm they match your intent.

4. Sensitivity analysis. Vary each input distribution and check whether the output changes materially. If it does, that input is driving the uncertainty and deserves more attention.

5. Compare to historical data. If you are simulating a metric you have historical data for, check that the simulated distribution is consistent with observed values.

A simulation is a model. It inherits all the risks of modeling: wrong assumptions, ignored correlations, inappropriate distributions. Validation is not optional — it is the difference between a useful analysis and a sophisticated error.