Monte Carlo Estimator

Last lecture, we built a language for randomness: events, random variables, expectation, variance, and covariance. Today, we use that vocabulary for the course’s first estimator: approximate an average with a sample average.

Most quantities in machine learning are expectations. The training loss we will minimize all semester is an average over a dataset far too large to touch every example at every step. The score of a game-playing agent is the expected outcome over the randomness of the game. Neither can be computed exactly, but both can be sampled. We need to quantify how well the average of a few random samples approximates the average over all possibilities. The error will be \(\sigma/\sqrt{n}\), where \(n\) is the number of samples and \(\sigma\) is the standard deviation of one sample.

Rolling with Advantage

Roll two fair six-sided dice and keep the larger value. Tabletop gamers call this rolling with advantage. It certainly feels better than rolling one die, but how much better, on average? Let \(X_1\) and \(X_2\) be the two rolls, the same two dice we met last lecture, so that a single game is the pair \(X = (X_1, X_2)\) and the value we keep is \(f(X) = \max(X_1, X_2)\). The number we are after is the long-run average value kept: \[ \mu = \mathbb{E}[\max(X_1, X_2)]. \]

For this tiny game, we can compute \(\mu\) exactly by enumerating the \(6 \times 6 = 36\) equally likely outcomes. The maximum is at most \(k\) exactly when both dice are at most \(k\), and the two rolls are independent, so those two conditions multiply: \[ \Pr(\max(X_1, X_2) \leq k) = \frac{k}{6} \cdot \frac{k}{6} = \frac{k^2}{36}. \] Subtracting consecutive values of that expression leaves the probability of landing on \(k\) exactly: \[ \Pr(\max(X_1, X_2) = k) = \frac{k^2 - (k-1)^2}{36} = \frac{2k - 1}{36}. \] Weighting each value by its probability gives the expectation we came for: \[ \mu = \sum_{k=1}^6 k \cdot \frac{2k-1}{36} = \frac{1 + 6 + 15 + 28 + 45 + 66}{36} = \frac{161}{36} \approx 4.47. \]

In the plot, advantage drags the distribution to the right: the flat gray distribution of a single die becomes a ramp that favors large values, and the mean climbs from \(3.5\) to about \(4.47\).

Enumeration worked because the game was tiny. Replace two dice with twenty, add rerolls and a rulebook of modifiers, and the number of outcomes explodes. Replace the game with “the error of a model on a random photo from the internet” and enumeration is not even meaningful. Even when we cannot enumerate the game, we can still play it.

The Sample Mean Estimator

Let \(X\) be a random variable (possibly a vector, like an image or a whole game), and let \(f\) score it with a real number. The quantity we want is the average score: \[ \mu = \mathbb{E}[f(X)]. \] We cannot evaluate that expectation, but we can draw \(n\) independent samples \(X^{(1)}, \ldots, X^{(n)}\), each distributed like \(X\), and report their sample mean. This rule is the Monte Carlo estimator: \[ \hat{\mu}_n = \frac{1}{n} \sum_{i=1}^n f(X^{(i)}). \] An estimator is any rule that turns samples into a guess for an unknown quantity; this one is named after the casino in Monaco. In the advantage roll, \(X^{(i)}\) is the \(i\)-th game, a fresh pair of dice, and \(f(X^{(i)})\) is the larger of that pair. The parenthesized superscript indexes samples, and it is the notation we will use for data points all semester.

We will check whether the estimator is correct on average and how its variance decreases as we collect samples.

Correct on Average

Linearity of expectation makes an average of \(n\) quantities with mean \(\mu\) itself equal to \(\mu\) on average, regardless of how the samples interact.

Claim: The Monte Carlo estimator is correct on average: \(\mathbb{E}[\hat{\mu}_n] = \mu\).

Proof of Claim

Each sample is distributed like \(X\), so \(\mathbb{E}[f(X^{(i)})] = \mu\) for every \(i\). By linearity of expectation, the constant \(\frac1n\) slides out and the expectation of the sum is the sum of the expectations: \[ \mathbb{E}[\hat{\mu}_n] = \frac{1}{n} \sum_{i=1}^n \mathbb{E}[f(X^{(i)})] = \frac{1}{n} \cdot n \mu = \mu. \] The proof did not use independence. Linearity of expectation holds no matter how the samples relate, so correlated samples can still give an unbiased estimate, though with higher variance.

An estimator whose expectation equals its target is called unbiased: its errors have no systematic direction. Unbiasedness alone is a low bar, though. The lazy rule “play one game and report it” is also unbiased, but nobody would trust it. What separates the sample mean from the lazy rule is the second question: variance.

The \(1/\sqrt{n}\) Rate

Variance is where independence finally matters. Let \(\sigma^2\) be the variance of a single sample’s score: \[ \sigma^2 = \textnormal{Var}(f(X)). \] Because the samples are independent, every covariance term in last lecture’s variance-of-a-sum formula is zero, so the \(n\) variances simply add. The factor \(\frac1n\) is squared by variance, so it divides the sum’s variance by \(n^2\): \[ \textnormal{Var}(\hat{\mu}_n) = \frac{\sigma^2}{n}. \] Problem 1 has you carry out this calculation in full, and then discover what happens to it when independence fails.

The error we should typically expect is the standard deviation, which is the square root of that variance: \[ \sqrt{\textnormal{Var}(\hat{\mu}_n)} = \frac{\sigma}{\sqrt{n}}. \] It shrinks to zero as \(n\) grows, which is the law of large numbers with a rate attached, and we call it the \(1/\sqrt{n}\) rate. The square root in the denominator is the bad news: halving the error takes four times the samples, and each extra decimal digit of accuracy takes a factor of one hundred more games. The numerator is the good news, because \(\sigma\) is the only feature of the problem that the formula mentions. How complicated \(f\) is, how many dice are on the table, how many dimensions \(X\) lives in: none of it enters except through the spread of a single game.

For the advantage roll, enumeration gives \(\sigma^2 = \frac{2555}{1296} \approx 1.97\), so \(\sigma \approx 1.40\), and the formula predicts an error near \(1.40/\sqrt{100} = 0.14\) after \(100\) games and near \(0.014\) after \(10{,}000\). The class demo plays the game at each \(n\) and checks that the measured error lands on this prediction.

Although one more digit costs a hundred times the samples, deterministic grids become much worse in high dimensions.

Any Dimension You Like

A \(1/\sqrt{n}\) rate may look slow until you see what it replaces. The deterministic alternative to sampling is a grid: to average a function over a square, evaluate it at every point of a \(20 \times 20\) lattice and average the \(400\) values. In two dimensions, that works fine. But the cost of a grid is \(m^d\) evaluations for \(m\) points per axis in \(d\) dimensions: at \(d = 10\) our modest grid already needs \(20^{10} \approx 10^{13}\) points, and the quantities we will meet later in the course live in millions of dimensions. This exponential blow-up is called the curse of dimensionality, and it kills the grid.

The Monte Carlo estimator does not: its error is \(\sigma/\sqrt{n}\), and that formula never mentions \(d\). A classic example is estimating \(\pi\) by throwing darts: let \(X\) be a uniformly random point in the unit square and \(f(X)\) the indicator that \(X\) lands inside the quarter circle, so that \(\mu = \mathbb{E}[f(X)] = \pi/4\). The exact same code estimates the volume of a ball in \(100\) dimensions, where no grid could ever go. (\(\sigma\) itself can grow with the problem; the rate \(1/\sqrt{n}\) does not.)

In the plot, both panels use \(400\) points to estimate the quarter circle’s area: a regular grid on the left and random darts on the right both land within a couple percent of \(\pi\).

In high dimensions, then, sampling is the only method still standing, so if we want a better answer, its error formula is what we have to attack.

Control Variate

Increasing \(n\) is expensive, but we can also reduce \(\sigma\), the spread of a single game. Suppose that alongside our target \(f(X)\) there is a surrogate quantity \(g(X)\) correlated with \(f(X)\) whose expectation \(\mathbb{E}[g(X)]\) is known exactly. For the advantage roll, the surrogate is already sitting on the table: take \(g(X) = X_1 + X_2\), the sum of the same two dice. Last lecture, linearity of expectation gave us \(\mathbb{E}[X_1 + X_2] = 7\) with no enumeration at all, and the sum is clearly correlated with the max, since a high sum usually means a high larger die.

The control variate adjustment subtracts the surrogate’s measured luck: \[ \hat{\mu}_c = \frac{1}{n} \sum_{i=1}^n \Big[ f(X^{(i)}) - c \big( g(X^{(i)}) - \mathbb{E}[g(X)] \big) \Big], \] where \(c\) is a constant we get to choose. (The subscript on \(\hat{\mu}_c\) names that coefficient rather than a sample count; we still draw \(n\) samples.) The correction term \(g(X^{(i)}) - \mathbb{E}[g(X)]\) measures which way the luck ran on game \(i\). If the dice summed above \(7\), the rolls ran lucky, the recorded max is probably an overestimate, and the adjustment pulls it back down by \(c\) times the measured luck.

The second in-class exercise checks that this random adjustment does not change the estimator’s expectation.

Claim: The adjusted estimator is unbiased for every choice of \(c\): \(\mathbb{E}[\hat{\mu}_c] = \mu\).

Proof of Claim

By linearity of expectation, we can take the expectation of each piece separately: \[ \mathbb{E}[\hat{\mu}_c] = \frac{1}{n} \sum_{i=1}^n \Big( \mathbb{E}[f(X^{(i)})] - c \big( \mathbb{E}[g(X^{(i)})] - \mathbb{E}[g(X)] \big) \Big) = \frac{1}{n} \sum_{i=1}^n \big( \mu - c \cdot 0 \big) = \mu, \] where each correction term vanishes because \(X^{(i)}\) is distributed like \(X\), so \(\mathbb{E}[g(X^{(i)})] = \mathbb{E}[g(X)]\). This is why the surrogate’s mean must be known: we subtract the true mean, not a guess, and the correction is centered exactly at zero.

The choice of \(c\) does not change the bias; it changes only the variance. Setting \(c = 0\) recovers the plain estimator, so the adjusted family can only do as well or better, provided we choose \(c\) well. The larger of two dice is roughly half their sum, so take \(c = \frac12\). Enumeration gives the variance of a single adjusted sample as \(\frac{665}{1296} \approx 0.51\), down from \(\sigma^2 \approx 1.97\). That is nearly a factor of four less variance, and since the error is the square root of the variance, it is half the error at every single \(n\): with \(n\) games the adjusted estimator reaches what the plain one needs about \(4n\) to match.

Problem 2 determines whether \(\frac12\) is optimal. You will show that the variance of \(\hat{\mu}_c\) is a parabola in \(c\), find its minimum in closed form, and show that the reduction is determined by the correlation \(\rho\) between \(f(X)\) and \(g(X)\), the covariance divided by the two standard deviations from last lecture. In Problem 1, positive \(\rho\) increased the poll’s variance. Here, the magnitude of \(\rho\) determines the variance reduction.

Variance Reduction Later in the Course

Today introduced a theme we return to for the rest of the semester: variance reduction. When we train neural networks, stochastic gradient descent will turn out to be a Monte Carlo estimator of the true gradient, with the batch playing the role of the samples. And in the final week, an agent learning from rewards will estimate its expected return by Monte Carlo and reduce the variance by subtracting a baseline, which is today’s adjustment down to the exact coefficient you derive on Problem 2.

Next lecture turns to linear algebra, the language in which data and models are written. When part of a measured quantity has a known average, subtracting that centered part can reduce the samples needed for a given error.