Monte Carlo Simulation
Advanced Statistical Methods
Solving Impossible Problems With Random Numbers
Monte Carlo methods use random sampling to approximate solutions to problems that are intractable analytically. The convergence rate of O(1/√n) is dimension-independent, making these methods powerful in high dimensions.
- Physics — Simulate particle interactions in nuclear reactor design and shielding calculations
- Finance — Price complex derivatives using risk-neutral Monte Carlo pricing models
- Engineering — Perform structural reliability analysis with probability of failure estimation
Monte Carlo methods harness randomness to tame the curse of dimensionality.
Importance Sampling
Rejection Sampling
import numpy as np
from scipy import stats
class MonteCarloSampler:
def __init__(self, seed=42):
self.rng = np.random.RandomState(seed)
def rejection_sampling(self, target_pdf, proposal_pdf, proposal_sampler, M, n_samples):
samples = []
attempts = 0
while len(samples) < n_samples:
x = proposal_sampler()
u = self.rng.uniform()
if u <= target_pdf(x) / (M * proposal_pdf(x)):
samples.append(x)
attempts += 1
return np.array(samples), len(samples) / attempts
def importance_sampling(self, target_pdf, proposal_pdf, proposal_sampler, n_samples, h=lambda x: x):
x = proposal_sampler(n_samples)
w = target_pdf(x) / proposal_pdf(x)
w_norm = w / np.sum(w)
mean = np.sum(w_norm * h(x))
var = np.sum(w_norm * (h(x) - mean)**2)
ess = np.sum(w)**2 / np.sum(w**2)
return mean, var, ess
def mcmc_gibbs(self, log_target, initial, n_samples, proposal_std=1.0):
n_vars = len(initial)
samples = np.zeros((n_samples, n_vars))
current = initial.copy()
accepted = np.zeros(n_vars)
for t in range(n_samples):
for i in range(n_vars):
candidate = current.copy()
candidate[i] += self.rng.normal(0, proposal_std)
log_ratio = log_target(candidate) - log_target(current)
if np.log(self.rng.uniform()) < log_ratio:
current = candidate
accepted[i] += 1
samples[t] = current
acceptance_rates = accepted / n_samples
return samples, acceptance_rates