Monte Carlo Basics

Monte Carlo methods estimate expectations and integrals by random sampling.

The simplest case is an expectation value

$$\mu = \mathbb E[f(X)].$$

If $X_1,\ldots,X_N$ are independent samples with the same distribution as $X$, then

$$\hat\mu_N = \frac1N \sum_{k=1}^N f(X_k)$$

is the basic Monte Carlo estimator. Its statistical error decreases like $N^{-1/2}$ under standard finite-variance assumptions.

Monte Carlo is not magic. It trades deterministic quadrature error for statistical uncertainty, and practical calculations must account for variance, bias, autocorrelation, sampling design, and reproducibility.

For how statistical error bars fit into broader numerical uncertainty budgets, see Error Estimates.

Expectations as Integrals

Many integrals can be written as expectation values. If $X$ has density $p(x)$, then

$$\mathbb E[f(X)] = \int f(x)p(x)\,dx.$$

Thus estimating an integral can be turned into estimating an average over random samples.

If the target integral is

$$I = \int_D g(x)\,dx,$$

and $D$ has finite volume $\operatorname{vol}(D)$, one may sample uniformly from $D$ and write

$$I = \operatorname{vol}(D)\, \mathbb E_{\mathrm{unif}}[g(X)].$$

The Monte Carlo estimate is then

$$\hat I_N = \operatorname{vol}(D) \frac1N \sum_{k=1}^N g(X_k).$$

This is often useful in high-dimensional problems where tensor-product quadrature grids become impossible.

Sample-Mean Estimator

Let

$$Y=f(X), \qquad \mu=\mathbb E[Y], \qquad \sigma_Y^2=\operatorname{Var}(Y).$$

For independent samples $Y_k=f(X_k)$, the sample mean

$$\bar Y_N = \frac1N \sum_{k=1}^N Y_k$$

is unbiased:

$$\mathbb E[\bar Y_N]=\mu.$$

Its variance is

$$\operatorname{Var}(\bar Y_N) = \frac{\sigma_Y^2}{N}.$$

The standard deviation of the estimator is therefore

$$\operatorname{SE}(\bar Y_N) = \frac{\sigma_Y}{\sqrt N}.$$

This standard error is the origin of the $1/\sqrt N$ Monte Carlo error scaling.

Estimated Error Bars

The true $\sigma_Y$ is usually unknown. From independent samples one estimates it with

$$s_N^2 = \frac{1}{N-1} \sum_{k=1}^N (Y_k-\bar Y_N)^2.$$

The estimated standard error is

$$\widehat{\operatorname{SE}} = \frac{s_N}{\sqrt N}.$$

For large $N$ and finite variance, the central limit theorem gives the approximate distribution

$$\frac{\bar Y_N-\mu}{\sigma_Y/\sqrt N} \approx \mathcal N(0,1).$$

This approximation justifies familiar error bars, but it can fail badly for heavy-tailed observables, strong correlations, or insufficient equilibration.

Error Scaling and Cost

The $N^{-1/2}$ rate is slow but dimension-independent in form. To reduce a Monte Carlo standard error by a factor of $10$, one usually needs about $100$ times as many independent samples.

Monte Carlo is attractive when:

• the dimension is high;
• sampling from the relevant distribution is easier than gridding the domain;
• modest stochastic accuracy is enough;
• deterministic quadrature is blocked by complex geometry or many degrees of freedom.

Monte Carlo is unattractive when:

• the integrand has huge variance;
• rare events dominate the answer;
• samples are strongly correlated;
• signs or phases cause severe cancellations;
• systematic bias dominates statistical error.

The rate $N^{-1/2}$ describes only statistical variance. It does not account for model error, discretization error, Markov-chain bias, floating-point error, or a poor estimator.

Importance Sampling

Importance sampling changes the sampling distribution to reduce variance or make sampling possible.

Suppose the goal is

$$I = \int f(x)p(x)\,dx,$$

but samples are drawn from a proposal density $q(x)$. If $q(x)>0$ wherever $f(x)p(x)$ contributes, then

$$I = \int f(x) \frac{p(x)}{q(x)} q(x)\,dx.$$

Thus

$$\hat I_N = \frac1N \sum_{k=1}^N f(X_k)w(X_k), \qquad w(x)=\frac{p(x)}{q(x)},$$

with $X_k\sim q$.

The weights correct the change of sampling distribution. Good importance sampling makes $f(x)w(x)$ less variable than direct sampling. Bad importance sampling can make variance enormous.

The support condition matters. If $q(x)=0$ in a region where $p(x)f(x)$ contributes, no finite weight can recover the missing region.

Self-Normalized Weights

Sometimes the target density is known only up to a normalization constant:

$$p(x) = \frac{\tilde p(x)}{Z}.$$

If samples come from $q$, a common estimator for

$$\mathbb E_p[f]$$

is the self-normalized importance estimate

$$\hat\mu_N = \frac{ \sum_{k=1}^N f(X_k)\tilde w(X_k) }{ \sum_{k=1}^N \tilde w(X_k) }, \qquad \tilde w(x)=\frac{\tilde p(x)}{q(x)}.$$

This estimator is generally biased at finite $N$, though often consistent under suitable conditions. Large variation in the weights is a warning sign that the proposal distribution has poor overlap with the target.

Markov Chains and Autocorrelation

Many useful distributions cannot be sampled independently. Markov-chain Monte Carlo constructs a correlated sequence

$$X_1,X_2,\ldots$$

whose stationary distribution is the desired target.

The same sample mean is often used, but the error bar changes. Positive autocorrelation reduces the effective number of independent samples:

$$N_{\mathrm{eff}} \approx \frac{N}{2\tau_{\mathrm{int}}},$$

where $\tau_{\mathrm{int}}$ is an integrated autocorrelation time. Correspondingly,

$$\operatorname{Var}(\bar Y_N) \approx \frac{\sigma_Y^2}{N_{\mathrm{eff}}}.$$

This is why Markov-chain calculations need burn-in checks, autocorrelation estimates, blocking or batching analyses, and independent-chain comparisons.

Bias and Variance

Monte Carlo error reports should separate statistical variance from bias.

Variance is random fluctuation from finite sampling. It can often be reduced by more samples, variance reduction, or a better proposal distribution.

Bias is a systematic shift. It can come from:

• sampling the wrong distribution;
• insufficient burn-in;
• a biased estimator;
• finite time-step or discretization error;
• using an approximate model;
• stopping an optimization based on noisy estimates.

Increasing $N$ reduces variance but does not automatically remove bias.

Quantum-Mechanics Uses

Quantum mechanics uses Monte Carlo ideas in several distinct ways.

• Repeated projective or POVM measurements produce samples from Born probabilities. Estimating an expectation value from many shots is a Monte Carlo estimate of a measurement distribution.
• Variational Monte Carlo rewrites a variational energy as an expectation over configurations sampled from $\lvert\Psi_\theta\rvert^2$.
• Path-integral and worldline Monte Carlo sample configurations or paths in imaginary time when a positive sampling weight is available.
• Stochastic wavefunction and quantum-trajectory methods sample individual trajectories whose ensemble reproduces an open-system master equation.

These uses share probability tools but differ in physics and algorithms. The reusable probability facts are here; the detailed computational methods belong to dedicated computational and many-body pages. For the existing physics preview, see Variational Monte Carlo Preview.

Sign and Phase Problems

Monte Carlo works best when the target can be interpreted as a nonnegative probability density. Quantum calculations often produce oscillatory phases or positive and negative weights.

If an integral is written as

$$I = \int A(x)s(x)p(x)\,dx,$$

where $p(x)\ge0$ is sampled and $s(x)$ is a sign or phase factor, then the estimator may involve an average sign or phase. When that average becomes very small, the relative error can grow exponentially with system size, time, or inverse temperature.

This sign or phase problem is not a minor implementation nuisance. It can be the central obstruction in real-time path integrals, frustrated spin systems, and fermionic many-body calculations.

The Sign Problem Preview owns the many-body reweighting identity, free-energy scaling of the average phase, basis dependence, sign-free structures, and complexity qualifications.

Common Mistakes

• Reporting a Monte Carlo number without an uncertainty estimate.
• Treating correlated Markov-chain samples as independent samples.
• Confusing stochastic error with systematic bias.
• Assuming $1/\sqrt N$ scaling solves poor overlap or rare-event sampling.
• Forgetting support conditions in importance sampling.
• Trusting self-normalized weights without checking weight degeneracy.
• Ignoring burn-in, autocorrelation time, and random-seed reproducibility.
• Treating a sign problem as merely a need for more patience.

Cross-Links

• Random Variables
• Probability Densities
• Expectation Values
• Variance and Covariance
• Conditional Probability
• Gaussian Distributions
• Fisher Information
• Error Estimates
• Born Rule
• Variational Monte Carlo Preview

References

• N. Metropolis and S. Ulam, “The Monte Carlo Method,” Journal of the American Statistical Association 44, 335–341, 1949.
• N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of State Calculations by Fast Computing Machines,” Journal of Chemical Physics 21, 1087–1092, 1953.
• J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods, Methuen, 1964.
• J. S. Liu, Monte Carlo Strategies in Scientific Computing, Springer, 2001.
• C. P. Robert and G. Casella, Monte Carlo Statistical Methods, 2nd ed., Springer, 2004.
• M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics, Oxford University Press, 1999.
• W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, “Quantum Monte Carlo simulations of solids,” Reviews of Modern Physics 73, 33–83, 2001.

Exercises

1. Show that the sample mean is unbiased for $\mu=\mathbb E[Y]$.

Solution

For independent identically distributed samples $Y_1,\ldots,Y_N$,

$$\bar Y_N = \frac1N \sum_{k=1}^N Y_k.$$

Linearity of expectation gives

$$\mathbb E[\bar Y_N] = \frac1N \sum_{k=1}^N \mathbb E[Y_k] = \frac1N \sum_{k=1}^N \mu = \mu.$$

2. Derive the variance of the sample mean for independent samples.

Solution

If each $Y_k$ has variance $\sigma_Y^2$ and the samples are independent, then

$$\operatorname{Var}(\bar Y_N) = \operatorname{Var} \left( \frac1N \sum_{k=1}^N Y_k \right).$$

Using independence,

$$\operatorname{Var}(\bar Y_N) = \frac{1}{N^2} \sum_{k=1}^N \operatorname{Var}(Y_k) = \frac{N\sigma_Y^2}{N^2} = \frac{\sigma_Y^2}{N}.$$

3. A simulation has independent samples with sample standard deviation $s_N=4$ and $N=10{,}000$. Estimate the standard error of the sample mean.

Solution

Use

$$\widehat{\operatorname{SE}} = \frac{s_N}{\sqrt N}.$$

Thus

$$\widehat{\operatorname{SE}} = \frac{4}{100} = 0.04.$$

4. Prove the basic importance-sampling identity.

Solution

Assume $q(x)>0$ wherever $f(x)p(x)$ contributes. Then

$$\begin{aligned} \int f(x)p(x)\,dx &= \int f(x) \frac{p(x)}{q(x)} q(x)\,dx\\ &= \mathbb E_q \left[ f(X)\frac{p(X)}{q(X)} \right]. \end{aligned}$$

The Monte Carlo estimator is the sample average of the random variable inside this expectation.

5. A two-outcome quantum measurement gives outcome $1$ with Born probability $p$. If $N$ independent shots are used to estimate $p$ by the observed frequency $\hat p$, what is $\operatorname{Var}(\hat p)$?

Solution

Each shot is a Bernoulli random variable with variance

$$p(1-p).$$

The frequency $\hat p$ is the sample mean of $N$ independent Bernoulli samples. Therefore

$$\operatorname{Var}(\hat p) = \frac{p(1-p)}{N}.$$

The standard error is

$$\sqrt{\frac{p(1-p)}{N}}.$$