Monte Carlo Methods

calendar Dec 12, 2024 · 2 min read · probability monte-carlo simulation numerical-methods  ·
Share on: linkedin copy

Interactive visualization of Monte Carlo methods for solving complex problems through random sampling.

Principle

Use random sampling to solve deterministic or stochastic problems.

Estimating π - Interactive Animation

Integration Concept

Monte Carlo integration works by randomly sampling points and determining the ratio that fall under the curve.

graph TD
    A[Start: Define function f(x) and bounds a, b] --> B[Generate N random points in rectangle]
    B --> C{Is point under curve?}
    C -->|Yes| D[Count as inside]
    C -->|No| E[Count as outside]
    D --> F{More samples?}
    E --> F
    F -->|Yes| B
    F -->|No| G[Calculate: Area × inside/total]
    G --> H[Result: Integral estimate]
    
    style A fill:#e1f5ff
    style H fill:#c8e6c9
    style C fill:#fff9c4

Algorithm:

  1. Define integration bounds $[a, b]$ and maximum function value $M$
  2. Generate random points $(x, y)$ where $x \in [a, b]$, $y \in [0, M]$
  3. Check if $y \leq f(x)$ (point is under curve)
  4. Estimate: $\int_a^b f(x)dx \approx (b-a) \times M \times \frac{\text{points under curve}}{\text{total points}}$

Convergence: Error decreases as $O(1/\sqrt{N})$ where $N$ is the number of samples.

Estimating π

 1import numpy as np
 2
 3def estimate_pi(n_samples=1000000):
 4    """Estimate π using Monte Carlo"""
 5    # Random points in [0,1] x [0,1]
 6    x = np.random.uniform(0, 1, n_samples)
 7    y = np.random.uniform(0, 1, n_samples)
 8    
 9    # Count points inside quarter circle
10    inside = (x**2 + y**2) <= 1
11    pi_estimate = 4 * np.mean(inside)
12    
13    return pi_estimate
14
15print(f"π ≈ {estimate_pi():.6f}")

Numerical Integration

$$ \int_a^b f(x) dx \approx \frac{b-a}{N} \sum_{i=1}^N f(x_i) $$
1def monte_carlo_integrate(f, a, b, n_samples=100000):
2    """Integrate f from a to b using Monte Carlo"""
3    x = np.random.uniform(a, b, n_samples)
4    return (b - a) * np.mean(f(x))
5
6# Example: integrate x^2 from 0 to 1
7result = monte_carlo_integrate(lambda x: x**2, 0, 1)
8print(f"Integral ≈ {result:.6f} (exact: 0.333333)")

Applications

  • Option pricing (finance)
  • Risk analysis
  • Bayesian inference (MCMC)
  • Physics simulations

Further Reading

Related Snippets