Monte Carlo Risk Model from Scratch

Point estimates lie. When someone says "the project will cost $2M," they're presenting one number as if the future is certain. A Monte Carlo model replaces that single number with a probability distribution — showing you not just the expected value but the full range of plausible outcomes and the likelihood of each.

This tutorial builds a project cost risk model using two Analysis Tool Box functions: CreateMetalogDistribution to fit a flexible distribution to your expert estimates, and SimulateNormallyDistributedOutcome for simpler, normally distributed components.

Prerequisites

Python 3.9+
analysistoolbox installed (pip install analysistoolbox)
Basic familiarity with probability distributions

What is a metalog distribution?

The metalog (metalogistic) distribution is a remarkably flexible continuous distribution that can take nearly any shape — symmetric, skewed, fat-tailed, bounded. It's fit from quantile estimates rather than from data, which makes it ideal for expert elicitation: you ask a subject matter expert "what's your 10th, 50th, and 90th percentile estimate?" and fit the metalog to those three points.

Note

The metalog was introduced by Thomas Keelin in 2016 specifically for risk analysis and decision analysis. It's now used by professional risk analysts at major consultancies as a replacement for triangle distributions and PERT distributions.

Step 1: Set up the scenario

Imagine you're estimating total cost for a software project with three cost components:

Development labor: expert says P10=$400K, P50=$600K, P90=$900K (right-skewed — cost overruns are more likely than savings)
Infrastructure: roughly normally distributed, mean=$150K, std=$30K
Licensing: roughly normally distributed, mean=$80K, std=$15K

Step 2: Fit the metalog for development labor

python

from analysistoolbox.simulations import CreateMetalogDistribution

labor_dist = CreateMetalogDistribution(
    quantile_values=[400_000, 600_000, 900_000],
    quantile_probabilities=[0.10, 0.50, 0.90],
    lower_bound=200_000,
    upper_bound=2_000_000,
    number_of_terms=3,
    plot_distribution=True,
)

CreateMetalogDistribution fits the metalog and plots the resulting PDF and CDF. Look at the plot — does the shape match your intuition about the cost distribution? If it looks implausibly flat or has a weird kink, try adjusting the quantile estimates or increasing number_of_terms.

Teaching Note

The lower_bound and upper_bound parameters constrain the distribution. For costs, a lower bound of $200K (minimum possible spend) and a wide upper bound prevent the metalog from assigning probability to physically impossible values like negative costs. Always set bounds when modeling quantities that have natural limits.

Step 3: Simulate 10,000 scenarios

python

import numpy as np

N = 10_000
np.random.seed(412)

# Sample from the metalog distribution using inverse CDF
uniform_samples = np.random.uniform(0, 1, N)
labor_samples = np.array([labor_dist.quantile(p) for p in uniform_samples])

Step 4: Simulate the normally distributed components

python

from analysistoolbox.simulations import SimulateNormallyDistributedOutcome

infra_samples = SimulateNormallyDistributedOutcome(
    mean=150_000,
    standard_deviation=30_000,
    number_of_simulations=N,
    random_seed=412,
    plot_simulation_results=False,
)

license_samples = SimulateNormallyDistributedOutcome(
    mean=80_000,
    standard_deviation=15_000,
    number_of_simulations=N,
    random_seed=413,
    plot_simulation_results=False,
)

Step 5: Aggregate and analyze

python

import pandas as pd
import matplotlib.pyplot as plt

total_cost = labor_samples + infra_samples + license_samples

results_df = pd.DataFrame({
    'labor': labor_samples,
    'infrastructure': infra_samples,
    'licensing': license_samples,
    'total': total_cost,
})

# Key statistics
print(f"P10 total cost: ${np.percentile(total_cost, 10):,.0f}")
print(f"P50 total cost: ${np.percentile(total_cost, 50):,.0f}")
print(f"P90 total cost: ${np.percentile(total_cost, 90):,.0f}")
print(f"Probability > $1.5M: {(total_cost > 1_500_000).mean():.1%}")

# Plot the distribution of total cost
plt.figure(figsize=(10, 5))
plt.hist(total_cost, bins=80, color='#2563EB', alpha=0.7, edgecolor='white')
plt.axvline(np.percentile(total_cost, 10), color='#16A34A', linestyle='--', label='P10')
plt.axvline(np.percentile(total_cost, 50), color='#D97706', linestyle='--', label='P50')
plt.axvline(np.percentile(total_cost, 90), color='#DC2626', linestyle='--', label='P90')
plt.xlabel('Total Project Cost ($)')
plt.ylabel('Frequency')
plt.title('Monte Carlo: Project Cost Distribution (10,000 simulations)')
plt.legend()
plt.tight_layout()
plt.show()

Interpreting the results

The output gives you answers no point estimate can provide:

"What's the budget we need to be 90% confident we won't exceed?" → P90 value
"What's the probability we exceed our $1.5M cap?" → directly readable from the simulation
"Which cost component contributes most to variance?" → compare the standard deviations of each component's samples

Teaching Note

The power of Monte Carlo is that it compounds uncertainty correctly. If development labor is uncertain AND infrastructure is uncertain, a simple sum of point estimates tells you nothing about tail risk. The simulation shows you the full joint distribution of total cost, including scenarios where both go wrong simultaneously. Always run at least 10,000 iterations — fewer iterations give unstable tail percentiles.

Full code

python

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from analysistoolbox.simulations import CreateMetalogDistribution, SimulateNormallyDistributedOutcome

N = 10_000
np.random.seed(412)

# Fit metalog for labor costs
labor_dist = CreateMetalogDistribution(
    quantile_values=[400_000, 600_000, 900_000],
    quantile_probabilities=[0.10, 0.50, 0.90],
    lower_bound=200_000,
    upper_bound=2_000_000,
    number_of_terms=3,
    plot_distribution=True,
)

# Sample labor via inverse CDF
uniform_samples = np.random.uniform(0, 1, N)
labor_samples = np.array([labor_dist.quantile(p) for p in uniform_samples])

# Simulate normally distributed components
infra_samples = SimulateNormallyDistributedOutcome(
    mean=150_000, standard_deviation=30_000, number_of_simulations=N,
    random_seed=412, plot_simulation_results=False,
)
license_samples = SimulateNormallyDistributedOutcome(
    mean=80_000, standard_deviation=15_000, number_of_simulations=N,
    random_seed=413, plot_simulation_results=False,
)

total_cost = labor_samples + infra_samples + license_samples

print(f"P10: ${np.percentile(total_cost, 10):,.0f}")
print(f"P50: ${np.percentile(total_cost, 50):,.0f}")
print(f"P90: ${np.percentile(total_cost, 90):,.0f}")
print(f"Prob > $1.5M: {(total_cost > 1_500_000).mean():.1%}")