Monte Carlo

Monte Carlo

Identity

Reference System Usage

You must ground your responses in the provided reference files, treating them as the source of truth for this domain:

• For Creation: Always consult `references/patterns.md`. This file dictates how things should be built. Ignore generic approaches if a specific pattern exists here.
• For Diagnosis: Always consult `references/sharp_edges.md`. This file lists the critical failures and "why" they happen. Use it to explain risks to the user.
• For Review: Always consult `references/validations.md`. This contains the strict rules and constraints. Use it to validate user inputs objectively.

Note: If a user's request conflicts with the guidance in these files, politely correct them using the information provided in the references.

Monte Carlo Simulation

Patterns

Basic Monte Carlo

Description

Fundamental Monte Carlo estimation

Example

import numpy as np from typing import Callable, Tuple, List from dataclasses import dataclass

@dataclass class MCResult: """Result of Monte Carlo simulation.""" mean: float std_error: float confidence_interval: Tuple[float, float] n_samples: int converged: bool

class MonteCarloEstimator: """ Basic Monte Carlo estimator with convergence checking. """

def __init__( self, target_std_error: float = 0.01, confidence_level: float = 0.95, max_samples: int = 1_000_000, batch_size: int = 10_000 ): self.target_std_error = target_std_error self.confidence_level = confidence_level self.max_samples = max_samples self.batch_size = batch_size

z-score for confidence interval

from scipy.stats import norm self.z = norm.ppf((1 + confidence_level) / 2)

def estimate( self, sampler: Callable[[], np.ndarray], estimand: Callable[[np.ndarray], float] ) -> MCResult: """ Estimate expected value via Monte Carlo.

Args: sampler: Function that returns random samples estimand: Function to compute from samples

Returns: MCResult with mean, uncertainty, and convergence info """ running_sum = 0.0 running_sum_sq = 0.0 n_total = 0

while n_total < self.max_samples:

Generate batch of samples

samples = sampler() values = np.array([estimand(s) for s in samples])

Update running statistics

batch_n = len(values) running_sum += np.sum(values) running_sum_sq += np.sum(values ** 2) n_total += batch_n

Compute current estimate

mean = running_sum / n_total variance = (running_sum_sq / n_total) - mean ** 2 std_error = np.sqrt(variance / n_total)

Check convergence

if std_error < self.target_std_error: return MCResult( mean=mean, std_error=std_error, confidence_interval=( mean - self.z std_error, mean + self.z std_error ), n_samples=n_total, converged=True )

Max samples reached without convergence

mean = running_sum / n_total variance = (running_sum_sq / n_total) - mean ** 2 std_error = np.sqrt(variance / n_total)

return MCResult( mean=mean, std_error=std_error, confidence_interval=( mean - self.z std_error, mean + self.z std_error ), n_samples=n_total, converged=False )

Example: Estimate pi

def sample_unit_square(): return np.random.uniform(-1, 1, size=(1000, 2))

def inside_circle(points): return np.mean(np.sum(points * 2, axis=1) <= 1) 4

mc = MonteCarloEstimator(target_std_error=0.001) result = mc.estimate(sample_unit_square, inside_circle) print(f"Pi estimate: {result.mean:.6f} +/- {result.std_error:.6f}")

Variance Reduction

Description

Techniques to reduce Monte Carlo variance

Example

import numpy as np from typing import Callable, Optional

class VarianceReduction: """Collection of variance reduction techniques."""

@staticmethod def antithetic_variates( sampler: Callable[[int], np.ndarray], estimand: Callable[[np.ndarray], float], n_samples: int ) -> tuple: """ Antithetic variates: use negative correlation.

For uniform U, also use 1-U. Reduces variance when estimand is monotonic. """

Generate primary samples

u = sampler(n_samples // 2)

Generate antithetic samples (1 - u for uniform)

u_anti = 1 - u

Evaluate on both

y = np.array([estimand(s) for s in u]) y_anti = np.array([estimand(s) for s in u_anti])

Average pairs

y_combined = (y + y_anti) / 2

mean = np.mean(y_combined) std_error = np.std(y_combined) / np.sqrt(len(y_combined))

return mean, std_error

@staticmethod def control_variates( samples: np.ndarray, y: np.ndarray, control: np.ndarray, control_mean: float ) -> tuple: """ Control variates: use known-mean correlated variable.

y_cv = y - c * (control - E[control]) c chosen to minimize variance. """

Optimal coefficient

cov_y_control = np.cov(y, control)[0, 1] var_control = np.var(control) c_optimal = cov_y_control / var_control

Adjusted estimator

y_cv = y - c_optimal * (control - control_mean)

mean = np.mean(y_cv) std_error = np.std(y_cv) / np.sqrt(len(y_cv))

Variance reduction ratio

var_ratio = np.var(y_cv) / np.var(y)

return mean, std_error, var_ratio

@staticmethod def importance_sampling( target_density: Callable[[np.ndarray], float], proposal_sampler: Callable[[int], np.ndarray], proposal_density: Callable[[np.ndarray], float], estimand: Callable[[np.ndarray], float], n_samples: int ) -> tuple: """ Importance sampling: sample from easier distribution.

E_p[f(x)] = E_q[f(x) * p(x) / q(x)]

Use when target is hard to sample but easy to evaluate. """

Sample from proposal

samples = proposal_sampler(n_samples)

Compute importance weights

weights = np.array([ target_density(s) / proposal_density(s) for s in samples ])

Weighted estimator

values = np.array([estimand(s) for s in samples]) weighted_values = values * weights

Self-normalized importance sampling (more stable)

mean = np.sum(weighted_values) / np.sum(weights)

Effective sample size (diagnostic)

ess = np.sum(weights) 2 / np.sum(weights 2)

return mean, ess

@staticmethod def stratified_sampling( sampler: Callable[[float, float, int], np.ndarray], estimand: Callable[[np.ndarray], float], n_strata: int, samples_per_stratum: int ) -> tuple: """ Stratified sampling: divide domain into strata.

Reduces variance by ensuring coverage of entire domain. """ stratum_means = []

for i in range(n_strata):

Stratum boundaries

low = i / n_strata high = (i + 1) / n_strata

Sample within stratum

samples = sampler(low, high, samples_per_stratum) values = np.array([estimand(s) for s in samples]) stratum_means.append(np.mean(values))

Overall mean (equal stratum weights)

mean = np.mean(stratum_means) std_error = np.std(stratum_means) / np.sqrt(n_strata)

return mean, std_error

Mcmc Sampling

Description

Markov Chain Monte Carlo for complex distributions

Example

import numpy as np from typing import Callable, Optional, List from dataclasses import dataclass

@dataclass class MCMCResult: """Result of MCMC sampling.""" samples: np.ndarray acceptance_rate: float effective_sample_size: float r_hat: Optional[float] = None # Gelman-Rubin diagnostic

class MetropolisHastings: """ Metropolis-Hastings MCMC sampler.

For sampling from distributions known up to normalizing constant. """

def __init__( self, log_target: Callable[[np.ndarray], float], proposal_std: float = 1.0, dim: int = 1 ): self.log_target = log_target self.proposal_std = proposal_std self.dim = dim

def sample( self, n_samples: int, initial: Optional[np.ndarray] = None, burn_in: int = 1000, thin: int = 1 ) -> MCMCResult: """ Generate samples using Metropolis-Hastings.

Args: n_samples: Number of samples to return initial: Starting point (default: zeros) burn_in: Samples to discard for equilibration thin: Keep every thin-th sample """ if initial is None: initial = np.zeros(self.dim)

current = initial.copy() current_log_p = self.log_target(current)

samples = [] n_accepted = 0 n_total = burn_in + n_samples * thin

for i in range(n_total):

Propose new state

proposal = current + np.random.normal( 0, self.proposal_std, size=self.dim ) proposal_log_p = self.log_target(proposal)

Accept/reject

log_alpha = proposal_log_p - current_log_p if np.log(np.random.uniform()) < log_alpha: current = proposal current_log_p = proposal_log_p n_accepted += 1

Store after burn-in, with thinning

if i >= burn_in and (i - burn_in) % thin == 0: samples.append(current.copy())

samples = np.array(samples)

return MCMCResult( samples=samples, acceptance_rate=n_accepted / n_total, effective_sample_size=self._compute_ess(samples) )

def _compute_ess(self, samples: np.ndarray) -> float: """Compute effective sample size from autocorrelation.""" n = len(samples) if samples.ndim > 1: samples = samples[:, 0] # Use first dimension

Compute autocorrelation

mean = np.mean(samples) var = np.var(samples) if var == 0: return n

autocorr = np.correlate(samples - mean, samples - mean, mode='full') autocorr = autocorr[n-1:] / (var * n)

Sum autocorrelations until negative

tau = 1.0 for k in range(1, n): if autocorr[k] < 0: break tau += 2 * autocorr[k]

return n / tau

class HamiltonianMonteCarlo: """ Hamiltonian Monte Carlo (HMC) sampler.

Uses gradient information for efficient exploration. """

def __init__( self, log_target: Callable[[np.ndarray], float], grad_log_target: Callable[[np.ndarray], np.ndarray], step_size: float = 0.1, n_leapfrog: int = 10, dim: int = 1 ): self.log_target = log_target self.grad_log_target = grad_log_target self.step_size = step_size self.n_leapfrog = n_leapfrog self.dim = dim

def sample( self, n_samples: int, initial: Optional[np.ndarray] = None, burn_in: int = 1000 ) -> MCMCResult: """Generate samples using HMC.""" if initial is None: initial = np.zeros(self.dim)

current_q = initial.copy() samples = [] n_accepted = 0

for i in range(burn_in + n_samples):

Sample momentum

current_p = np.random.normal(size=self.dim)

Leapfrog integration

q = current_q.copy() p = current_p.copy()

Half step for momentum

p += 0.5 self.step_size self.grad_log_target(q)

Full steps

for _ in range(self.n_leapfrog - 1): q += self.step_size p p += self.step_size self.grad_log_target(q)

Final half step

q += self.step_size p p += 0.5 self.step_size * self.grad_log_target(q)

Negate momentum for reversibility

p = -p

Compute Hamiltonian

current_H = -self.log_target(current_q) + 0.5 np.sum(current_p 2) proposed_H = -self.log_target(q) + 0.5 np.sum(p ** 2)

Accept/reject

if np.log(np.random.uniform()) < current_H - proposed_H: current_q = q n_accepted += 1

if i >= burn_in: samples.append(current_q.copy())

return MCMCResult( samples=np.array(samples), acceptance_rate=n_accepted / (burn_in + n_samples), effective_sample_size=len(samples) # HMC typically has high ESS )

Quasi Monte Carlo

Description

Low-discrepancy sequences for faster convergence

Example

import numpy as np from typing import Generator

class QuasiMonteCarlo: """ Quasi-Monte Carlo using low-discrepancy sequences.

Converges as O(1/N) vs O(1/sqrt(N)) for standard MC. """

@staticmethod def halton_sequence(dim: int, n: int, skip: int = 0) -> np.ndarray: """ Generate Halton sequence.

Low-discrepancy sequence using prime bases. """ def halton_single(index: int, base: int) -> float: result = 0.0 f = 1.0 / base i = index while i > 0: result += f * (i % base) i //= base f /= base return result

First dim primes

primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37][:dim]

points = np.zeros((n, dim)) for i in range(n): for d in range(dim): points[i, d] = halton_single(i + skip + 1, primes[d])

return points

@staticmethod def sobol_sequence(dim: int, n: int) -> np.ndarray: """ Generate Sobol sequence.

Better uniformity than Halton for high dimensions. """ from scipy.stats import qmc sampler = qmc.Sobol(d=dim, scramble=True) return sampler.random(n)

@staticmethod def latin_hypercube(dim: int, n: int) -> np.ndarray: """ Latin Hypercube Sampling.

Ensures each dimension has uniform marginal coverage. """ from scipy.stats import qmc sampler = qmc.LatinHypercube(d=dim) return sampler.random(n)

Comparison: QMC vs MC for integration

def integrate_comparison(f, dim, n_points):

Standard Monte Carlo

mc_samples = np.random.uniform(size=(n_points, dim)) mc_estimate = np.mean([f(s) for s in mc_samples])

Quasi Monte Carlo (Sobol)

qmc_samples = QuasiMonteCarlo.sobol_sequence(dim, n_points) qmc_estimate = np.mean([f(s) for s in qmc_samples])

return mc_estimate, qmc_estimate

Sensitivity Analysis

Description

Global sensitivity analysis with Monte Carlo

Example

import numpy as np from typing import Callable, Dict, List from dataclasses import dataclass

@dataclass class SobolIndices: """Sobol sensitivity indices.""" first_order: Dict[str, float] # S_i total_order: Dict[str, float] # ST_i second_order: Dict[tuple, float] # S_ij

class SensitivityAnalysis: """ Global sensitivity analysis using Sobol indices.

Quantifies how much each input contributes to output variance. """

def __init__( self, model: Callable[[np.ndarray], float], param_bounds: Dict[str, tuple], n_samples: int = 10000 ): self.model = model self.param_names = list(param_bounds.keys()) self.bounds = np.array(list(param_bounds.values())) self.n_samples = n_samples self.dim = len(self.param_names)

def compute_sobol_indices(self) -> SobolIndices: """ Compute first-order and total-order Sobol indices.

Uses Saltelli's extension of Sobol sequence. """

Generate two independent sample matrices

from scipy.stats import qmc sampler = qmc.Sobol(d=self.dim, scramble=True)

A = sampler.random(self.n_samples) B = sampler.random(self.n_samples)

Scale to bounds

A = self._scale_samples(A) B = self._scale_samples(B)

Evaluate model on base matrices

f_A = np.array([self.model(a) for a in A]) f_B = np.array([self.model(b) for b in B])

Compute indices for each parameter

first_order = {} total_order = {}

for i, name in enumerate(self.param_names):

Matrix with column i from B, rest from A

AB_i = A.copy() AB_i[:, i] = B[:, i] f_AB_i = np.array([self.model(ab) for ab in AB_i])

Matrix with column i from A, rest from B

BA_i = B.copy() BA_i[:, i] = A[:, i] f_BA_i = np.array([self.model(ba) for ba in BA_i])

First-order index

var_total = np.var(np.concatenate([f_A, f_B])) first_order[name] = np.mean(f_B * (f_AB_i - f_A)) / var_total

Total-order index

total_order[name] = 0.5 np.mean((f_A - f_AB_i) * 2) / var_total

return SobolIndices( first_order=first_order, total_order=total_order, second_order={} # Computed separately if needed )

def _scale_samples(self, samples: np.ndarray) -> np.ndarray: """Scale [0,1] samples to parameter bounds.""" return samples * (self.bounds[:, 1] - self.bounds[:, 0]) + self.bounds[:, 0]

def morris_screening(self, n_trajectories: int = 10) -> Dict[str, tuple]: """ Morris method for parameter screening.

Cheaper than Sobol, identifies important vs unimportant parameters. """ from scipy.stats import qmc

delta = 1.0 / (2 * (n_trajectories - 1)) elementary_effects = {name: [] for name in self.param_names}

for _ in range(n_trajectories):

Random starting point

x = np.random.uniform(size=self.dim) * (1 - delta)

Random permutation of parameters

order = np.random.permutation(self.dim)

x_scaled = self._scale_samples(x.reshape(1, -1))[0] f_current = self.model(x_scaled)

for i in order: x_new = x.copy() x_new[i] += delta

x_new_scaled = self._scale_samples(x_new.reshape(1, -1))[0] f_new = self.model(x_new_scaled)

Elementary effect

ee = (f_new - f_current) / delta elementary_effects[self.param_names[i]].append(ee)

x = x_new f_current = f_new

Compute mean and std of elementary effects

return { name: (np.mean(np.abs(ee)), np.std(ee)) for name, ee in elementary_effects.items() }

Anti-Patterns

---

Pattern

Fixed sample size without convergence check

Problem

May waste samples or stop before convergence

Solution

Use adaptive sampling with target error threshold

---

Pattern

Single random seed for all runs

Problem

Results not reproducible, can't parallelize properly

Solution

Use separate streams: np.random.SeedSequence for spawning

---

Pattern

Ignoring autocorrelation in MCMC

Problem

Effective samples much less than nominal

Solution

Compute ESS, use thinning or better sampler (HMC)

---

Pattern

Uniform sampling for peaked functions

Problem

Most samples contribute nothing to estimate

Solution

Use importance sampling with proposal near peak

---

Pattern

No burn-in for MCMC

Problem

Samples biased by initial state

Solution

Discard burn-in samples, check convergence diagnostics

---

Pattern

High-dimensional Halton sequence

Problem

Correlations appear in high dimensions

Solution

Use Sobol with scrambling for d > 10

Monte Carlo - Sharp Edges

Insufficient Samples for Rare Events

Id

not-enough-samples

Severity

critical

Summary

Standard MC fails for low-probability events

Symptoms

• Estimated probability is 0 but event can happen
• Huge variance in rare event estimates
• Different runs give wildly different results

Why

Standard MC samples uniformly. For probability p, need ~1/p samples to see one event. For p = 10^-6, need millions of samples.

Worse: variance of rare event estimate is p(1-p)/N ~ p/N Relative error = 1/sqrt(Np) - huge for small p.

Critical for: safety analysis, risk assessment, extreme events.

Gotcha

Estimate probability of rare failure

def simulate_failure(): return some_complex_model() < threshold # Returns True rarely

n_samples = 10000 failures = sum(simulate_failure() for _ in range(n_samples)) p_failure = failures / n_samples # Often 0!

"Zero failures in 10000 runs" doesn't mean failure is impossible

Just means p < ~1/10000

Solution

1. Importance sampling - sample near failure region

def importance_sample_failure():

Sample from distribution concentrated near failure

biased_sample = sample_near_threshold() is_failure = evaluate(biased_sample)

Correct for bias with likelihood ratio

weight = true_density(biased_sample) / biased_density(biased_sample) return is_failure * weight

2. Subset simulation - chain of conditional probabilities

class SubsetSimulation: def __init__(self, p0: float = 0.1): self.p0 = p0 # Conditional probability per level

def estimate(self, model, n_samples: int) -> float: samples = [sample_prior() for _ in range(n_samples)] responses = [model(s) for s in samples]

levels = [] current_threshold = np.percentile(responses, (1-self.p0)*100) levels.append(current_threshold)

while current_threshold > target_threshold:

MCMC to sample conditional distribution

(samples with response > current_threshold)

samples = self.mcmc_conditional(samples, responses, current_threshold) responses = [model(s) for s in samples] current_threshold = np.percentile(responses, (1-self.p0)*100) levels.append(current_threshold)

Probability = p0^(n_levels) * fraction above final threshold

n_above = sum(1 for r in responses if r > target_threshold) return self.p0 * len(levels) n_above / n_samples

3. Cross-entropy method for rare events

Inadequate Random Number Generator

Id

poor-random-numbers

Severity

high

Summary

Weak RNG introduces bias or correlation

Symptoms

• Results depend on initial seed in unexpected ways
• Parallel runs produce correlated results
• Patterns appear in 'random' samples

Why

Not all RNGs are created equal. Linear congruential generators (LCG): short period, correlations. Many default RNGs fail statistical tests.

For MC simulation:

• Need long period (> samples you'll ever use)
• Need independence across dimensions
• Need reproducibility for debugging
• Parallel streams must not overlap

Using time-based seeds in parallel = disaster.

Gotcha

import random random.seed() # Uses system time

Parallel processes started same second = same seed = same results!

Or worse: overlapping sequences give correlated estimates

Old-style seeding

np.random.seed(42)

Global state - other code can change it!

Solution

1. Use modern RNGs with proper seeding

from numpy.random import Generator, PCG64

Per-thread generators from seed sequence

from numpy.random import SeedSequence

ss = SeedSequence(12345) child_seeds = ss.spawn(n_processes) generators = [Generator(PCG64(s)) for s in child_seeds]

Each process uses its own generator

def worker(rng: Generator): samples = rng.uniform(size=(1000, dim))

Guaranteed independent from other workers

2. For reproducibility

rng = Generator(PCG64(seed=42)) # Local state, not global

3. Check RNG quality

Use Dieharder or TestU01 for serious applications

4. For cryptographic/security: use secrets module

MCMC Chain Not Converged to Target Distribution

Id

mcmc-not-converged

Severity

high

Summary

Samples from wrong distribution due to insufficient burn-in

Symptoms

• Results depend strongly on starting point
• Multiple chains give different answers
• Posterior mean drifts over long runs

Why

MCMC needs time to "forget" its starting point. Chain must reach stationary distribution. No guaranteed time - depends on problem.

Signs of non-convergence:

• Trace plots show trend
• Different chains disagree
• Gelman-Rubin R-hat > 1.1

Can't prove convergence, only detect non-convergence.

Gotcha

Run chain without checking convergence

samples = mcmc.sample(n_samples=10000, burn_in=100) # Burn-in too short!

First samples still influenced by starting point

Posterior mean is biased

Or: chain stuck in local mode

Never explores full posterior

Solution

1. Run multiple chains from dispersed starting points

def run_with_diagnostics(n_chains: int = 4, n_samples: int = 10000): chains = [] for i in range(n_chains): initial = sample_overdispersed_prior() chain = mcmc.sample(n_samples, initial=initial) chains.append(chain)

Gelman-Rubin diagnostic

r_hat = compute_r_hat(chains) if r_hat > 1.1: raise ValueError(f"Chains not converged: R-hat = {r_hat}")

return np.concatenate(chains)

2. Check trace plots

def visual_diagnostics(chain): import matplotlib.pyplot as plt fig, axes = plt.subplots(2, 1) axes[0].plot(chain) # Trace plot - should be stationary axes[1].hist(chain, bins=50) # Should match expected posterior

3. Compute effective sample size

ess = compute_ess(chain) if ess < 100: print(f"Warning: ESS only {ess}, need more samples")

4. Use better samplers (HMC, NUTS) for faster mixing

Monte Carlo Fails in High Dimensions

Id

curse-of-dimensionality

Severity

high

Summary

Convergence rate doesn't improve but constant gets worse

Symptoms

• Need exponentially more samples as dimension grows
• Integration error grows with dimension
• QMC loses advantage over MC

Why

MC convergence: O(1/sqrt(N)), independent of dimension. Sounds good, but the constant grows exponentially!

In d dimensions:

• Volume of hypercube: 1
• Volume of inscribed hypersphere: approaches 0
• Most samples are in "corners", far from center

For integration: most samples contribute nothing. For MCMC: harder to explore space efficiently.

Gotcha

High-dimensional integration

def integrate_gaussian(dim): samples = np.random.uniform(-3, 3, size=(100000, dim)) values = np.exp(-np.sum(samples * 2, axis=1) / 2) return np.mean(values) (6 ** dim)

d=2: works well

d=10: huge variance

d=100: all samples give ~0, estimate is garbage

Solution

1. Dimension reduction if possible

Find low-dimensional structure in problem

2. Use QMC with dimensionality-aware sampling

from scipy.stats import qmc

Sobol better than Halton for high-d

sampler = qmc.Sobol(d=dim, scramble=True) samples = sampler.random(n)

3. Importance sampling in important dimensions

Identify which dimensions matter most

4. Sequential Monte Carlo

Build up distribution through tempering

5. Tensor decomposition methods

Exploit structure to reduce effective dimension

6. Active subspaces

Find low-dim subspace where function varies

Biased Monte Carlo Estimator

Id

estimator-bias

Severity

medium

Summary

Estimator systematically over/underestimates true value

Symptoms

• Error doesn't decrease with more samples
• Confidence interval excludes true value
• Systematic offset in estimates

Why

Not all MC estimators are unbiased. Common sources of bias:

1. Self-normalized importance sampling E[sum(w*f)/sum(w)] != E[f] in general

2. Ratio estimators E[X/Y] != E[X]/E[Y]

3. Nonlinear functions of expectations E[g(X)] != g(E[X]) unless g linear

4. Truncated or censored samples

Gotcha

Biased ratio estimator

def estimate_ratio(): x_samples = simulate_x(n=1000) y_samples = simulate_y(n=1000) return np.mean(x_samples) / np.mean(y_samples) # Biased!

Self-normalized importance sampling

def snis_estimate(samples, weights, f_values): return np.sum(weights * f_values) / np.sum(weights) # Biased!

Solution

1. Use unbiased estimators when possible

def unbiased_ratio_estimator(n):

Generate paired samples

pairs = [(simulate_x(), simulate_y()) for _ in range(n)]

Use f(x)/g(y) with proper weights

return np.mean([x / y for x, y in pairs if y != 0])

2. For self-normalized IS: use more samples

Bias is O(1/n), variance is O(1/n)

Bias becomes negligible with enough samples

3. Bias correction

def jackknife_bias_correction(estimates): n = len(estimates) theta_all = np.mean(estimates) theta_i = [(np.sum(estimates) - e) / (n - 1) for e in estimates] bias = (n - 1) * (np.mean(theta_i) - theta_all) return theta_all - bias

4. Bootstrap confidence intervals

Account for bias in interval construction

Importance Sampling Weights Have Infinite Variance

Id

variance-explosion

Severity

medium

Summary

Heavy-tailed weights cause unstable estimates

Symptoms

• Occasional samples dominate entire estimate
• ESS is tiny compared to sample size
• Estimates jump wildly between runs

Why

Importance sampling weight: w(x) = p(x)/q(x)

If q(x) has lighter tails than p(x):

• w(x) can be huge in tails
• Var(w) may be infinite
• CLT doesn't apply, convergence is slow

One bad sample can dominate thousands of good ones.

Gotcha

Sampling Normal(0,1) with Uniform(-5,5) proposal

def bad_importance_sampling(): samples = np.random.uniform(-5, 5, size=10000)

Normal density at x=4: very small

Uniform density at x=4: 0.1

Weight = small/0.1 = still small

But in tails: Normal decays faster than any polynomial

Uniform is constant

Weights are well-behaved here...

Worse: sampling heavy-tailed with light-tailed proposal

def worse_importance_sampling():

Sample Cauchy with Normal proposal

samples = np.random.normal(size=10000)

Cauchy tails: O(1/x^2)

Normal tails: O(exp(-x^2))

Ratio = exp(x^2) / x^2 -> infinity!

A few extreme samples = entire estimate

Solution

1. Use proposal with heavier tails than target

If target is Normal, proposal can be t-distribution

2. Truncate weights

def truncate_weights(weights, max_weight=None): if max_weight is None: max_weight = 10 * np.median(weights) return np.minimum(weights, max_weight)

3. Self-normalized estimator (reduces variance)

def self_normalized(samples, weights, f): return np.sum(weights * f(samples)) / np.sum(weights)

4. Monitor effective sample size

def effective_sample_size(weights): w_normalized = weights / np.sum(weights) return 1 / np.sum(w_normalized ** 2)

If ESS << n, proposal is poor

5. Adaptive importance sampling

Update proposal based on observed samples

Monte Carlo - Validations

Using Global Random State

Id

global-random-state

Severity

warning

Type

regex

Pattern

• np\.random\.seed\(|random\.seed\(
• np\.random\.(rand|randn|uniform|normal)\(

Message

Global random state is not thread-safe and hard to reproduce. Use Generator.

Fix Action

Use rng = np.random.Generator(np.random.PCG64(seed)) with rng.uniform()

Applies To

• */.py

Fixed Sample Size Without Convergence Check

Id

fixed-sample-size

Severity

info

Type

regex

Pattern

• for\s+_\s+in\s+range\(\d+\):\s\n.sample
• n_samples\s=\s\d+(?!.*converge|error|std)

Message

Fixed sample size may waste computation or stop before convergence.

Fix Action

Use adaptive sampling with target error: while std_error > target: sample more

Applies To

• */.py

MCMC Without Burn-in

Id

no-burn-in-mcmc

Severity

warning

Type

regex

Pattern

• mcmc\.sample\(.\)(?!.burn|warm)
• metropolis.\(.\)(?!.*burn)

Message

MCMC chains need burn-in to reach equilibrium.

Fix Action

Add burn_in parameter: mcmc.sample(n_samples, burn_in=1000)

Applies To

• */.py

MCMC Without Effective Sample Size Check

Id

no-ess-check

Severity

info

Type

regex

Pattern

• samples\s=.mcmc(?!.*ess|effective)

Message

MCMC effective sample size may be much less than nominal.

Fix Action

Compute ESS: if ess < 100: warn or run longer

Applies To

• */.py

Importance Sampling Without Weight Diagnostics

Id

importance-no-weight-check

Severity

warning

Type

regex

Pattern

• weight.=.density.\/.density(?!.*ess|max|var)
• importance_weight(?!.*check|valid|ess)

Message

Importance weights can have huge variance. Check ESS and max weights.

Fix Action

Check: ess = 1/sum(w_norm**2); if ess < n/10: proposal is poor

Applies To

• */.py

Halton Sequence in High Dimensions

Id

halton-high-dim

Severity

info

Type

regex

Pattern

• halton.dim\s=\s[1-9]\d+|halton.d\s=\s[1-9]\d+

Message

Halton sequences show correlation patterns in d > 10. Use Sobol instead.

Fix Action

For high dimensions: sampler = qmc.Sobol(d=dim, scramble=True)

Applies To

• */.py

Parallel MC Without Proper Seeding

Id

no-seed-parallel

Severity

warning

Type

regex

Pattern

• Pool\(.\)\.map.random(?!.*SeedSequence|spawn)
• parallel.random(?!.seed.*sequence)

Message

Parallel random streams may overlap without proper seeding.

Fix Action

Use SeedSequence.spawn() for independent parallel streams

Applies To

• */.py

MCMC With Single Chain Only

Id

single-chain-mcmc

Severity

info

Type

regex

Pattern

• mcmc\.sample\(.\)(?!.chain|n_chain|multiple)

Message

Single MCMC chain can't detect convergence issues. Run multiple chains.

Fix Action

Run 4+ chains from dispersed starts, compute R-hat diagnostic

Applies To

• */.py

MC Estimate Without Error Bound

Id

no-variance-estimate

Severity

info

Type

regex

Pattern

• np\.mean\(.sample.\)(?!.*std|var|error|ci)

Message

Monte Carlo estimates should include uncertainty quantification.

Fix Action

Compute std error: stderr = np.std(samples) / np.sqrt(n)

Applies To

• */.py

QMC Sequence Without Scrambling

Id

qmc-without-scramble

Severity

info

Type

regex

Pattern

• Sobol\(.scramble\s=\s*False
• Halton\((?!.*scramble)

Message

Unscrambled QMC sequences can show bias. Enable scrambling.

Fix Action

Use scramble=True: qmc.Sobol(d=dim, scramble=True)

Applies To

• */.py