Monte Carlo Methods & Variance Reduction

Monte Carlo methods estimate a quantity by sampling rather than exhaustive summation or deterministic quadrature. They are especially valuable when dimension makes a regular grid impossible: doubling resolution in each of $d$ dimensions multiplies grid work by $2^d$, whereas ordinary Monte Carlo's root-mean-square error remains proportional to $N^{-1/2}$ under finite variance. That rate is slow, but it does not deteriorate directly with dimension.

Estimation begins with an expectation #

Suppose the target is

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

For independent draws $X_i\sim p$,

$$\hat I_N=\frac1N\sum_{i=1}^N f(X_i)$$

is unbiased, with variance $\mathop{\mathrm{Var}}_p[f(X)]/N$. An estimated Monte Carlo standard error (MCSE) is

$$\widehat{\mathop{\mathrm{MCSE}}} =\frac{s_f}{\sqrt N},$$

where $s_f$ is the sample standard deviation of the evaluated quantity. Reporting only $\hat I_N$ hides numerical uncertainty from finite simulation.

To reduce independent-sampling error by a factor of ten requires roughly one hundred times as many draws. Better algorithms therefore try to reduce variance per unit cost, not merely increase $N$.

Variance reduction uses structure #

Antithetic and stratified sampling #

Antithetic draws couple samples so their errors tend to cancel; for a monotone function under a uniform input, pairing $U$ with $1-U$ is a common example. Stratification divides the input space into regions and samples each, preventing random omission of important strata. Both require a structure that actually correlates with $f$.

Control variates #

If $g(X)$ has known expectation $\mu_g$, estimate

$$I=\mathbb E[f(X)-\beta(g(X)-\mu_g)].$$

The optimal scalar $\beta=\mathop{\mathrm{Cov}}(f,g)/\mathop{\mathrm{Var}}(g)$ reduces variance when $g$ tracks $f$. Estimating $\beta$ from the same data needs care in finite samples.

Importance sampling #

When direct draws from $p$ are hard but a proposal $q$ is easy,

$$I=\mathbb E_q\left[f(X)\frac{p(X)}{q(X)}\right].$$

The proposal must cover every region where $f(x)p(x)$ matters. A few extreme weights create a formally correct but practically unstable estimator. Inspect normalized weights, tail behavior, and effective sample size; a large raw sample does not compensate for missing support.

MCMC trades independence for accessible transitions #

Markov-chain Monte Carlo constructs a transition kernel whose invariant distribution is the target $\pi$. After initialization and finite transitions, states are correlated and may still reflect the start. For a scalar estimand, the variance of the sample mean is approximately

$$\mathop{\mathrm{Var}}(\bar f)\approx \frac{\mathop{\mathrm{Var}}_\pi(f)}{N_{\mathrm{eff}}},$$

where $N_{\mathrm{eff}}$ accounts for autocorrelation. ESS and MCSE belong to a particular parameter or function; there is no single chain-wide sample quality number.

Run multiple dispersed chains. Rank-normalized split-$\hat R$ compares between- and within-chain behavior and can reveal nonstationarity or scale disagreement; bulk and tail ESS quantify different regions. Values near one for $\hat R$ are necessary evidence, not proof that all modes were found. A shared initialization in one mode can make every chain agree on the wrong region.

Warmup adapts and moves chains toward typical regions; discarding a finite prefix cannot guarantee convergence. Thinning usually throws away information and should be used for storage constraints, not as a cure for autocorrelation. A better transition kernel or longer run is the remedy.

Gibbs sampling uses full conditional distributions #

For a joint target $\pi(x_1,\ldots,x_d)$, a systematic Gibbs sweep draws

$$\begin{aligned} x_1'&\sim\pi(x_1\mid x_2,\ldots,x_d),\\ x_2'&\sim\pi(x_2\mid x_1',x_3,\ldots,x_d),\\ &\ \vdots\\ x_d'&\sim\pi(x_d\mid x_1',\ldots,x_{d-1}'). \end{aligned}$$

Each update uses the newest available coordinates. To see why one coordinate update preserves the joint target, suppose $(X,Y)\sim\pi$. Keep $Y$ and redraw $X'\sim\pi(x\mid Y)$. Then

$$\Pr(X'\in A,Y\in B) =\int_B\pi(A\mid y)\pi_Y(dy) =\pi(A\times B).$$

Composing invariant coordinate kernels yields an invariant sweep. This explains correctness at stationarity; it does not establish irreducibility or fast mixing from an arbitrary start.

Block Gibbs samples several strongly coupled variables jointly. Random-scan Gibbs selects coordinates randomly; systematic scan applies a fixed order. Both can preserve the target under appropriate conditions, but they can have different reversibility and convergence rates.

Worked bivariate-normal Gibbs sampler #

For a standard bivariate normal with correlation $\rho$,

$$\begin{bmatrix}X\\Y\end{bmatrix} \sim\mathcal N\left( \begin{bmatrix}0\\0\end{bmatrix}, \begin{bmatrix}1&\rho\\\rho&1\end{bmatrix} \right),$$

the full conditionals are

$$X\mid Y=y\sim\mathcal N(\rho y,1-\rho^2), \qquad Y\mid X=x\sim\mathcal N(\rho x,1-\rho^2).$$

When $|\rho|$ is small, alternating draws move freely. As $|\rho|\to1$, the conditional variance shrinks and the chain advances in small correlated steps along a narrow diagonal. The target is unchanged, but ESS per sweep collapses. Joint rotation, blocking, or a different sampler addresses geometry that more iterations only brute-force.

Gibbs in Boltzmann machines #

An RBM's bipartite graph makes hidden units independent given the visible layer and visible units independent given the hidden layer. One block Gibbs step can therefore update an entire layer in parallel. Contrastive divergence begins near data and uses only a short chain, so its negative sample is intentionally not a converged independent draw from the model. That bias is an optimization approximation and must not be described as exact maximum likelihood.

A reproducible Monte Carlo report #

Record the target, estimand, sampler and transition parameters, initialization, warmup/adaptation, saved draws, independent chains, random seeds, ESS, MCSE, $\hat R$, trace/autocorrelation checks, and any divergences or mode-specific behavior. For importance sampling, add proposal support and weight diagnostics. For Gibbs, add scan order, blocks, and conditional parameterization.

Simulation produces an estimate plus simulation error. A method is not validated by producing many numbers; it is validated by showing that its draws explore the intended target and that the remaining Monte Carlo uncertainty is small enough for the decision.