One-line summary: Apply Monte Carlo integration, MCMC sampling, and statistical mechanics simulations using numpy, scipy, and numba for physics research.
When to Use This Skill
When computing high-dimensional integrals analytically intractable (partition functions, path integrals)
When sampling from complex probability distributions (Boltzmann, posterior distributions)
When simulating statistical mechanical systems (Ising model, lattice gauge theory)
When estimating uncertainties via bootstrap or error propagation
When implementing Markov Chain Monte Carlo (Metropolis-Hastings, Gibbs sampling)
When studying phase transitions and critical phenomena
Trigger keywords: Monte Carlo integration, MCMC, Metropolis-Hastings, Ising model, importance sampling, statistical mechanics, partition function, random sampling
Background & Key Concepts
Monte Carlo Integration
Monte Carlo integration estimates an integral by random sampling:
where samples $\mathbf{x}_i \sim p(\mathbf{x})$. The statistical error scales as $1/\sqrt{N}$ regardless of dimensionality — Monte Carlo's key advantage over quadrature rules in high dimensions.
Importance Sampling
For efficiency, sample from a proposal distribution $q(\mathbf{x})$ and reweight:
$$
\langle f \rangle_p = \frac{\langle f \cdot w \rangle_q}{\langle w \rangle_q}, \quad w(\mathbf{x}) = \frac{p(\mathbf{x})}{q(\mathbf{x})}
$$
Choose $q$ to be large where $|f \cdot p|$ is large.
Markov Chain Monte Carlo
MCMC constructs a Markov chain whose stationary distribution is the target $p(\mathbf{x})$. The Metropolis-Hastings acceptance probability:
$$
H = -J \sum_{\langle i,j \rangle} s_i s_j - h \sum_i s_i
$$
where $s_i \in {-1, +1}$ are spins, $J$ is the coupling constant, and $h$ is external field. The critical temperature is $T_c = 2J / k_B \ln(1 + \sqrt{2}) \approx 2.269 J/k_B$.
Example 1: Quantum Harmonic Oscillator Ground State Energy
# =============================================# Path integral Monte Carlo: harmonic oscillator# =============================================import numpy as np
import matplotlib.pyplot as plt
defpath_integral_harmonic(N_beads=50, N_steps=100_000, beta=5.0, omega=1.0,
hbar=1.0, m=1.0, seed=42):
"""
Estimate ground state energy of harmonic oscillator via PIMC.
E_0 = hbar*omega/2 = 0.5 (in natural units).
"""
rng = np.random.default_rng(seed)
tau = beta / N_beads
dtau = tau
path = rng.normal(0, 1/np.sqrt(m*omega), N_beads)
energies = []
for step inrange(N_steps):
# Random bead update
k = rng.integers(N_beads)
delta = rng.normal(0, 0.3)
x_new = path[k] + delta
# Periodic boundary: bead k-1 and k+1
k_m = (k - 1) % N_beads
k_p = (k + 1) % N_beads
# Action difference
dS = (m / (2 * dtau)) * ((x_new - path[k_m])**2 - (path[k] - path[k_m])**2 +
(path[k_p] - x_new)**2 - (path[k_p] - path[k])**2)
dS += dtau * 0.5 * m * omega**2 * (x_new**2 - path[k]**2)
if rng.uniform() < np.exp(-dS):
path[k] = x_new
if step > N_steps // 4:
# Virial estimator
E = 0.5 / beta + 0.5 * m * omega**2 * np.mean(path**2)
energies.append(E)
energies = np.array(energies)
E_mean = energies.mean()
E_err = energies.std() / np.sqrt(len(energies))
return E_mean, E_err
E, dE = path_integral_harmonic()
print(f"E_0 (PIMC) = {E:.4f} ± {dE:.4f}")
print(f"E_0 (exact) = 0.5000")
Interpreting these results: The PIMC estimate should converge to 0.5 ħω. Increase N_beads and N_steps to reduce statistical and systematic errors.
Example 2: Bootstrap Error Estimation for Experimental Data
# =============================================# Bootstrap resampling for uncertainty quantification# =============================================import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(42)
# Simulate experimental measurements with noise
N_exp = 200
true_mean, true_std = 3.14, 0.5
data = rng.normal(true_mean, true_std, N_exp)
defbootstrap_ci(data, statistic_fn, n_boot=5000, ci=95, seed=42):
"""
Bootstrap confidence interval for any statistic.
"""
rng = np.random.default_rng(seed)
boot_stats = np.empty(n_boot)
n = len(data)
for b inrange(n_boot):
sample = rng.choice(data, size=n, replace=True)
boot_stats[b] = statistic_fn(sample)
lo = np.percentile(boot_stats, (100 - ci) / 2)
hi = np.percentile(boot_stats, 100 - (100 - ci) / 2)
return boot_stats, lo, hi
# Estimate mean and its 95% CI
boot_means, lo, hi = bootstrap_ci(data, np.mean)
print(f"Sample mean: {data.mean():.4f}")
print(f"Bootstrap 95% CI: [{lo:.4f}, {hi:.4f}]")
print(f"Analytical 95% CI: [{data.mean() - 1.96*data.std()/np.sqrt(N_exp):.4f}, "f"{data.mean() + 1.96*data.std()/np.sqrt(N_exp):.4f}]")
# Estimate skewness and its CI
boot_skew, lo_sk, hi_sk = bootstrap_ci(data, stats.skew)
print(f"\nSkewness estimate: {stats.skew(data):.4f}")
print(f"Bootstrap 95% CI: [{lo_sk:.4f}, {hi_sk:.4f}]")
Interpreting these results: Bootstrap CIs are distribution-free — use for any statistic (correlation, skewness, custom physics observable) where analytical formulas are unavailable.
Last updated: 2026-03-17 | Maintainer: @xjtulycIssues: GitHub Issues
