Monte Carlo Simulation Guide

The Monte Carlo engine in stochlab provides high-performance parallel simulation with advanced features for large-scale experiments.

Quick Start

from stochlab.models import MarkovChain
from stochlab.mc import MonteCarloEngine
import numpy as np

# Create a process
P = np.array([[0.7, 0.3], [0.4, 0.6]])
mc = MarkovChain.from_transition_matrix(["A", "B"], P)

# Create engine and simulate
engine = MonteCarloEngine(mc)
result = engine.simulate(n_paths=10000, T=100, seed=42)

print(f"Simulated {len(result.paths)} paths")

Core Concepts

What is Monte Carlo Simulation?

Monte Carlo simulation estimates quantities by:

Running a stochastic process many times (generating sample paths)
Computing a statistic on each path
Averaging the results

Example: Estimate the probability that a Markov chain reaches state B by time 100.

def reaches_b(path):
    return 1.0 if "B" in path.states else 0.0

stats = engine.estimate(
    estimator_fn=reaches_b,
    n_paths=10000,
    T=100,
    parallel=True
)

print(f"P(reach B) = {stats.mean:.4f} ± {stats.stderr:.4f}")
print(f"95% CI: {stats.confidence_interval}")

Features

1. Parallel Simulation

Run simulations in parallel to leverage multiple CPU cores:

# Sequential (uses 1 core)
result = engine.simulate(n_paths=10000, T=100, parallel=False)

# Parallel (uses all cores)
result = engine.simulate(n_paths=10000, T=100, parallel=True)

# Parallel with specific number of workers
result = engine.simulate(n_paths=10000, T=100, parallel=True, n_jobs=4)

Performance: Expect 6-8x speedup on an 8-core machine.

When to use parallel:

n_paths > 100 (overhead is amortized)
Long time horizons (T is large)
Complex processes (each path takes time)

When NOT to use parallel:

Small simulations (n_paths < 50)
Very fast processes (overhead dominates)
Memory-constrained environments

2. Memory-Efficient Modes

Choose how much data to store based on your needs:

Paths Mode (Default)

Store complete path trajectories:

result = engine.simulate(n_paths=1000, T=100, mode="paths")

# Access full paths
for path in result.paths:
    print(path.times)   # [0, 1, 2, ..., 100]
    print(path.states)  # ['A', 'B', 'A', ...]

Use when: You need complete trajectory information for visualization or detailed analysis.

Memory: ~10-20 KB per path (depends on path length).

Values Mode

Store only final state values:

result = engine.simulate(n_paths=1000, T=100, mode="values", parallel=True)

# Access final values
final_states = [path.states[-1] for path in result.paths]

Use when: You only care about endpoints (e.g., steady-state analysis).

Memory: ~8 bytes per path (90% reduction).

Stats Mode

Store only aggregated statistics:

result = engine.simulate(n_paths=1000, T=100, mode="stats", parallel=True)

# Statistics are in metadata
stats = result.metadata["statistics"]

Use when: Estimating expectations and don’t need individual paths.

Memory: ~100 bytes total (99% reduction).

3. Reproducible Seeds

Monte Carlo simulations are reproducible when you specify a seed:

# Run 1
result1 = engine.simulate(n_paths=1000, T=100, seed=42)

# Run 2 (identical results)
result2 = engine.simulate(n_paths=1000, T=100, seed=42)

# Verify
assert all(np.array_equal(p1.states, p2.states) 
           for p1, p2 in zip(result1.paths, result2.paths))

How it works: Uses NumPy’s SeedSequence to generate statistically independent seeds for each path while maintaining reproducibility.

4. Batch Size Control

The engine automatically computes optimal batch sizes, but you can override:

# Automatic (recommended)
result = engine.simulate(n_paths=10000, T=100, parallel=True)

# Manual control
result = engine.simulate(
    n_paths=10000,
    T=100,
    parallel=True,
    batch_size=100  # Paths per batch
)

Batch size trade-offs:

Smaller batches (50-100): Better load balancing, more overhead
Larger batches (500-1000): Lower overhead, worse load balancing
Optimal (100-200): Balance between the two

The automatic calculation targets 4-10 batches per worker, which works well for most cases.

5. Progress Tracking

Show progress bars for long-running simulations (requires tqdm):

result = engine.simulate(
    n_paths=100000,
    T=100,
    parallel=True,
    show_progress=True  # Requires: pip install tqdm
)

If tqdm is not installed, the simulation continues without a progress bar.

Advanced Usage

Estimating Expectations

The estimate() method computes expectations with confidence intervals:

def estimator(path):
    """Compute some quantity on the path."""
    return path.states[-1] == "B"  # Indicator function

stats = engine.estimate(
    estimator_fn=estimator,
    n_paths=10000,
    T=100,
    parallel=True,
    confidence_level=0.95
)

print(f"Estimate: {stats.mean:.4f}")
print(f"Std Error: {stats.stderr:.4f}")
print(f"95% CI: [{stats.confidence_interval[0]:.4f}, {stats.confidence_interval[1]:.4f}]")

Output:

Estimate: 0.5431
Std Error: 0.0050
95% CI: [0.5333, 0.5529]

Complex Estimators

You can compute any statistic:

# Average state occupancy
def avg_b_occupancy(path):
    return sum(s == "B" for s in path.states) / len(path.states)

# First passage time
def first_passage_time(path):
    try:
        return np.where(path.states == "B")[0][0]
    except IndexError:
        return float('inf')  # Never reached

# Path length until absorption
def absorption_time(path):
    return len(path) - 1

Working with Results

result = engine.simulate(n_paths=1000, T=100, parallel=True)

# Basic info
print(f"Number of paths: {len(result.paths)}")
print(f"Parallel: {result.metadata['parallel']}")
print(f"Batches used: {result.metadata.get('n_batches', 'N/A')}")

# Convert to DataFrame
df = result.to_dataframe()
print(df.head())
#    path_id  t  time state
# 0        0  0     0     A
# 1        0  1     1     B
# 2        0  2     2     B
# ...

# Analyze state distribution at specific time
dist = result.state_distribution(t=50)
print(dist)  # {'A': 0.42, 'B': 0.58}

Performance Tips

1. Use Parallel for Large Simulations

# Sequential: ~10 seconds
result = engine.simulate(n_paths=100000, T=100, parallel=False)

# Parallel (8 cores): ~1.5 seconds (6-7x speedup)
result = engine.simulate(n_paths=100000, T=100, parallel=True)

2. Choose the Right Mode

# If you need full paths
result = engine.simulate(n_paths=10000, T=1000, mode="paths")
# Memory: ~170 MB

# If you only need final values
result = engine.simulate(n_paths=10000, T=1000, mode="values")
# Memory: ~80 KB (2000x reduction!)

# If you only need statistics
result = engine.simulate(n_paths=10000, T=1000, mode="stats")
# Memory: ~1 KB (170,000x reduction!)

3. Don’t Over-Parallelize

# BAD: Too few paths for parallelism
result = engine.simulate(n_paths=10, T=100, parallel=True)
# Overhead dominates, slower than sequential!

# GOOD: Enough paths to amortize overhead
result = engine.simulate(n_paths=10000, T=100, parallel=True)
# Overhead is negligible, get full speedup

Rule of thumb: Use parallel=True when n_paths > 100.

4. Batch Size Matters

The default automatic batch sizing works well, but for very large simulations:

# For maximum throughput
result = engine.simulate(
    n_paths=1_000_000,
    T=100,
    parallel=True,
    batch_size=1000  # Larger batches = lower overhead
)

# For better progress tracking
result = engine.simulate(
    n_paths=1_000_000,
    T=100,
    parallel=True,
    batch_size=100,  # More frequent updates
    show_progress=True
)

Common Patterns

Pattern 1: Probability Estimation

# Estimate P(event occurs)
def event_occurs(path):
    return 1.0 if some_condition(path) else 0.0

stats = engine.estimate(event_occurs, n_paths=10000, T=100, parallel=True)
print(f"P(event) = {stats.mean:.4f} ± {stats.stderr:.4f}")

Pattern 2: Expected Value

# Estimate E[X_T]
def final_value(path):
    return float(path.states[-1])

stats = engine.estimate(final_value, n_paths=10000, T=100, parallel=True)
print(f"E[X_100] = {stats.mean:.4f}")

Pattern 3: Hitting Time Distribution

# Estimate distribution of hitting times
result = engine.simulate(n_paths=10000, T=1000, parallel=True)

hitting_times = []
for path in result.paths:
    hit_indices = np.where(path.states == "B")[0]
    if len(hit_indices) > 0:
        hitting_times.append(hit_indices[0])

print(f"Mean hitting time: {np.mean(hitting_times):.2f}")
print(f"Median hitting time: {np.median(hitting_times):.2f}")

Pattern 4: Comparing Processes

# Compare two different processes
engine1 = MonteCarloEngine(process1)
engine2 = MonteCarloEngine(process2)

result1 = engine1.simulate(n_paths=10000, T=100, seed=42, parallel=True)
result2 = engine2.simulate(n_paths=10000, T=100, seed=42, parallel=True)

# Compare distributions
dist1 = result1.state_distribution(t=100)
dist2 = result2.state_distribution(t=100)

Troubleshooting

Issue: Parallel simulation crashes on Windows

Symptom:

RuntimeError: An attempt has been made to start a new process...

Solution: Wrap your code in if __name__ == "__main__":

if __name__ == "__main__":
    engine = MonteCarloEngine(process)
    result = engine.simulate(n_paths=10000, T=100, parallel=True)

This is required on Windows due to how multiprocessing spawns processes.

Issue: Out of memory

Symptom: System runs out of RAM with large simulations.

Solution 1: Use values or stats mode:

result = engine.simulate(
    n_paths=1_000_000,
    T=1000,
    mode="values",  # Much lower memory
    parallel=True
)

Solution 2: Process in chunks:

results = []
for i in range(10):
    chunk = engine.simulate(
        n_paths=100000,
        T=1000,
        seed=42 + i,
        parallel=True
    )
    # Process chunk immediately
    results.append(process_chunk(chunk))

Issue: Slower with parallelization

Symptom: parallel=True is slower than parallel=False.

Causes:

Too few paths: Overhead dominates
Very fast processes: Each path takes <1ms
Small time horizon: T < 10

Solution: Use parallel=False for small simulations, or increase n_paths.

Issue: Results not reproducible

Symptom: Same seed gives different results.

Cause: Not specifying a seed, or process has internal randomness not controlled by RNG.

Solution: Always specify seed parameter:

result = engine.simulate(n_paths=10000, T=100, seed=42)  # Reproducible

API Reference

MonteCarloEngine

MonteCarloEngine(process: StochasticProcess)

Create a Monte Carlo engine for a stochastic process.

Parameters:

process: Any StochasticProcess (MarkovChain, RandomWalk, MM1Queue, etc.)

simulate()

simulate(
    n_paths: int,
    T: int,
    x0: State | None = None,
    *,
    parallel: bool = False,
    n_jobs: int = -1,
    batch_size: int | None = None,
    mode: str = "paths",
    seed: int | None = None,
    show_progress: bool = False,
    **kwargs
) -> SimulationResult

Run Monte Carlo simulation.

Parameters:

n_paths: Number of paths to simulate (must be > 0)
T: Time horizon for each path
x0: Initial state (optional, uses process default if None)
parallel: Enable parallel execution (default: False)
n_jobs: Number of workers, -1 = all CPUs (default: -1)
batch_size: Paths per batch (default: auto-computed)
mode: Return mode - “paths”, “values”, or “stats” (default: “paths”)
seed: Random seed for reproducibility (default: random)
show_progress: Show progress bar (default: False, requires tqdm)
**kwargs: Additional arguments passed to process.sample_path()

Returns: SimulationResult with paths and metadata

estimate()

estimate(
    estimator_fn: Callable[[Path], float],
    n_paths: int,
    T: int,
    x0: State | None = None,
    *,
    parallel: bool = False,
    confidence_level: float = 0.95,
    seed: int | None = None,
    **kwargs
) -> MCStatistics

Estimate E[f(X)] using Monte Carlo.

Parameters:

estimator_fn: Function that takes a Path and returns a scalar
n_paths: Number of Monte Carlo samples
T: Time horizon
x0: Initial state (optional)
parallel: Use parallel execution (default: False)
confidence_level: Confidence level for CI (default: 0.95)
seed: Random seed (optional)
**kwargs: Additional arguments for sample_path()

Returns: MCStatistics with mean, std, confidence interval, etc.

MCStatistics

Result from estimate() containing:

Attributes:

n_paths: Number of paths simulated
mean: Estimated expectation
std: Standard deviation of estimator
stderr: Standard error (std / sqrt(n))
confidence_interval: Tuple (lower, upper)
confidence_level: Confidence level (e.g., 0.95)
variance_reduction_factor: VRF if variance reduction used
metadata: Additional information

ReturnMode

Enum for simulation return modes:

ReturnMode.PATHS: Full Path objects
ReturnMode.VALUES: Final values only
ReturnMode.STATS: Aggregated statistics only

Summary

The Monte Carlo engine provides:

✓ High performance: 6-8x speedup with parallelization
✓ Memory efficiency: 90-99% reduction with values/stats modes
✓ Reproducibility: SeedSequence for quality randomness
✓ Ease of use: Simple API with sensible defaults
✓ Flexibility: Works with any stochastic process

For most users, this is all you need:

engine = MonteCarloEngine(process)
result = engine.simulate(n_paths=10000, T=100, parallel=True, seed=42)

Happy simulating! 🎲

Monte Carlo Simulation Guide

Quick Start

Core Concepts

What is Monte Carlo Simulation?

Features

1. Parallel Simulation

2. Memory-Efficient Modes

Paths Mode (Default)

Values Mode

Stats Mode

3. Reproducible Seeds

4. Batch Size Control

5. Progress Tracking

Advanced Usage

Estimating Expectations

Complex Estimators

Working with Results

Performance Tips

1. Use Parallel for Large Simulations

2. Choose the Right Mode

3. Don’t Over-Parallelize

4. Batch Size Matters

Common Patterns

Pattern 1: Probability Estimation

Pattern 2: Expected Value

Pattern 3: Hitting Time Distribution

Pattern 4: Comparing Processes

Troubleshooting

Issue: Parallel simulation crashes on Windows

Issue: Out of memory

Issue: Slower with parallelization

Issue: Results not reproducible

API Reference

MonteCarloEngine

simulate()

estimate()

MCStatistics

ReturnMode

Further Reading

Summary