Implicit Control Variate
8 min read
We want the average of a large population, but getting every value would be too expensive. Instead, we sample a few values and use them to estimate the average. We have two choices: whether to sample with or without replacement, and whether to use a fixed or random number of samples.
Sampling without replacement reduces the variance, especially when the sample approaches the whole population. Fixing the number of samples does something more surprising: it can reduce the expected squared error by a large constant factor. We’ll explore how.
Sampling Choices
A population holds \(n>1\) real values \(x_1,\dots,x_n\), with mean and variance: \[ \begin{gathered} \bar{x}=\frac1n\sum_{i=1}^n x_i,\\ \sigma^2=\frac1n\sum_{i=1}^n(x_i-\bar{x})^2. \end{gathered} \]
Let the budget \(m\) be an integer between \(1\) and \(n\). Every estimator below divides the sampled sum by \(m\), counting repeated draws separately. For a fixed count this is the empirical average; for a random count, \(m\) is the expected count and the denominator stays fixed, including when the sample is empty. All four estimators are unbiased, so their expected squared errors equal their variances:
| Sampling | Fixed Count: Exactly \(m\) | Random Count: \(m\) on Average |
|---|---|---|
| With Replacement | \(\displaystyle\frac{\sigma^2}{m}\) | \(\displaystyle\frac{\sigma^2+\bar{x}^{\,2}}{m}\) |
| Without Replacement | \(\displaystyle\frac{n-m}{n-1}\frac{\sigma^2}{m}\) | \(\displaystyle\left(1-\frac mn\right)\frac{\sigma^2+\bar{x}^{\,2}}{m}\) |
In the random-count column, sampling with replacement means taking a Poisson number of independent uniform draws. Sampling without replacement means independently including each population value with probability \(m/n\). The count is then binomial. The fixed-count column uses independent uniform draws with replacement or a uniformly chosen subset of exactly \(m\) distinct values.
Fixing the count removes the extra \(\bar{x}^2\) from the with-replacement variance. The plot shows the resulting gap for a population whose values share a common offset. All four curves use the same population and the same actual or expected budget.
The vertical scale is logarithmic above \(0.001\) and linear near zero, so the plot includes the exact zero errors of both without-replacement designs at \(m=n\). To see why the squared mean disappears, we can first subtract a constant by hand, then make that adjustment cancel by fixing the count.
Random Count
Draw a count \(K\) from a Poisson distribution with mean \(m\), independently of the sampled values. A Poisson count has equal mean and variance: \[ \mathbb{E}[K]=m, \qquad \mathrm{Var}(K)=m. \]
Each draw \(Y_j\) independently equals \(x_i\) with probability \(1/n\). It therefore has mean \(\bar{x}\) and variance \(\sigma^2\). Taking \(K\) such draws, with replacement, gives the estimate: \[ \hat{x}^{\mathrm{random}}=\frac1m\sum_{j=1}^K Y_j. \]
An empty sum is zero. If every value equals a nonzero \(a\), the estimate is \(aK/m\): it fluctuates whenever the sample size does, even though there is no variation in the population.
Given \(K\), the expected sampled sum is \(K\bar{x}\), so the estimate is unbiased: \[ \mathbb{E}[\hat{x}^{\mathrm{random}}] =\frac{\mathbb{E}[K]\bar{x}}m =\bar{x}. \]
Claim: The random-count estimator has variance: \[ \mathrm{Var}(\hat{x}^{\mathrm{random}}) =\frac{\sigma^2+\bar{x}^{\,2}}m. \]
Proof of the variance
Once the count is given, the draws are independent, so their variances add: \[\begin{align} \mathrm{Var}(\hat{x}^{\mathrm{random}}\mid K) &=\frac1{m^2}\sum_{j=1}^K\mathrm{Var}(Y_j\mid K)\\ &=\frac1{m^2}\sum_{j=1}^K\sigma^2\\ &=\frac{K\sigma^2}{m^2}. \end{align}\]
The law of total variance separates the average variance at a given count from the variance of the conditional mean: \[\begin{align} \mathrm{Var}(\hat{x}^{\mathrm{random}}) &=\mathbb{E}\!\left[\mathrm{Var}(\hat{x}^{\mathrm{random}}\mid K)\right] +\mathrm{Var}\!\left(\mathbb{E}[\hat{x}^{\mathrm{random}}\mid K]\right)\\ &=\mathbb{E}\!\left[\frac{K\sigma^2}{m^2}\right] +\mathrm{Var}\!\left(\frac{K\bar{x}}m\right)\\ &=\frac{\sigma^2}{m^2}\mathbb{E}[K] +\frac{\bar{x}^{\,2}}{m^2}\mathrm{Var}(K)\\ &=\frac{\sigma^2}{m^2}m+\frac{\bar{x}^{\,2}}{m^2}m\\ &=\frac{\sigma^2+\bar{x}^{\,2}}m. \end{align}\]
The second line substitutes the two conditional moments. The third pulls out the constants, and the fourth uses the Poisson mean and variance.
The first term comes from which values we draw. The second comes from how many we draw: the conditional mean moves with \(K\) whenever the population mean is nonzero. When the values cluster around a large common offset, that offset can dominate the variance without increasing the population’s spread.
Constant Adjustment
Subtracting a baseline lets us remove the common offset. Pick a constant \(c\), estimate the average of the shifted values \(Y_j-c\), and add \(c\) back at the end: \[ \hat{x}^{\mathrm{random}}_c =c+\frac1m\sum_{j=1}^K(Y_j-c). \]
For every fixed \(c\), the adjustment preserves unbiasedness: \[ \mathbb{E}[\hat{x}^{\mathrm{random}}_c] =c+\frac{\mathbb{E}[K]}m(\bar{x}-c) =\bar{x}. \]
Claim: The adjusted estimator has variance: \[ \mathrm{Var}(\hat{x}^{\mathrm{random}}_c) =\frac{\sigma^2+(\bar{x}-c)^2}m. \]
Proof
Shifting every draw by \(c\) changes its mean to \(\bar{x}-c\) and leaves its variance at \(\sigma^2\). Adding the constant back has no effect on variance. The preceding variance calculation, now with mean \(\bar{x}-c\), gives: \[\begin{align} \mathrm{Var}(\hat{x}^{\mathrm{random}}_c) &=\mathrm{Var}\!\left(\frac1m\sum_{j=1}^K(Y_j-c)\right)\\ &=\frac{\sigma^2}{m^2}\mathbb{E}[K] +\frac{(\bar{x}-c)^2}{m^2}\mathrm{Var}(K)\\ &=\frac{\sigma^2+(\bar{x}-c)^2}m. \end{align}\]
The last line again substitutes the Poisson mean and variance.
The variance is a parabola in \(c\), minimized when \(c\) is the unknown average itself: \[ \mathrm{Var}(\hat{x}^{\mathrm{random}}_{\bar{x}}) =\frac{\sigma^2}m. \]
Tuning \(c\) to \(\bar{x}\) removes the squared mean from the variance. Unfortunately, the optimal baseline is the very average we are trying to estimate. We cannot choose it exactly without already knowing the answer. Computing it from the whole population would defeat the reason we sampled in the first place.
In practice we could use an independent pilot sample to estimate the baseline, or update it from past samples. Either approach gives an approximation to the best adjustment. Subtracting a running reward baseline is also used in REINFORCE, a policy-gradient method in reinforcement learning [Williams, ML ’92].
Collecting the \(c\) terms shows what the baseline does: \[\begin{align} \hat{x}^{\mathrm{random}}_c &=c+\frac1m\sum_{j=1}^K Y_j-\frac1m\sum_{j=1}^K c\\ &=c+\hat{x}^{\mathrm{random}}-\frac{cK}m\\ &=\hat{x}^{\mathrm{random}}-\frac cm(K-m). \end{align}\]
The first line distributes the sum, the second counts the \(K\) copies of \(c\), and the third collects the two \(c\) terms. The baseline multiplies \(K-m\), the deviation of the sample size from its mean. It is therefore a control variate on the sample size: a correlated, mean-zero correction that we can subtract without changing the estimator’s expectation. A design with a fixed sample size makes this correction identically zero.
Fixing the Count
The beautiful part is that we can get exactly the same variance without knowing \(\bar{x}\). We only need to fix the number of samples. Draw exactly \(m\) values uniformly with replacement and take their empirical average: \[ \hat{x}^{\mathrm{fixed}}=\frac1m\sum_{j=1}^m Y_j. \]
Only the count has changed; values can still repeat. Now subtracting any baseline and adding it back cancels on every sample: \[\begin{align} \hat{x}^{\mathrm{fixed}}_c &=c+\frac1m\sum_{j=1}^m(Y_j-c)\\ &=c+\hat{x}^{\mathrm{fixed}}-\frac1m\sum_{j=1}^m c\\ &=c+\hat{x}^{\mathrm{fixed}}-\frac{mc}m\\ &=\hat{x}^{\mathrm{fixed}}. \end{align}\]
The second line splits off the empirical average, the third counts the \(m\) copies of \(c\), and the last cancels the adjustment. The corrected estimate is identically the plain one, for every \(c\). That includes the optimal \(c=\bar{x}\), which we no longer need to know.
Each draw has mean \(\bar{x}\), so the fixed-count estimate is unbiased: \[ \mathbb{E}[\hat{x}^{\mathrm{fixed}}] =\frac1m\sum_{j=1}^m\mathbb{E}[Y_j] =\bar{x}. \]
Claim: The fixed-count estimator has variance: \[ \mathrm{Var}(\hat{x}^{\mathrm{fixed}})=\frac{\sigma^2}m. \]
Proof of the variance
The draws are independent and each has variance \(\sigma^2\): \[\begin{align} \mathrm{Var}(\hat{x}^{\mathrm{fixed}}) &=\frac1{m^2}\sum_{j=1}^m\mathrm{Var}(Y_j)\\ &=\frac1{m^2}\sum_{j=1}^m\sigma^2\\ &=\frac{m\sigma^2}{m^2}\\ &=\frac{\sigma^2}m. \end{align}\]
The first line uses independence, the second substitutes the single-draw variance, and the last two count the terms and cancel a factor of \(m\).
Fixing the count achieves exactly the variance of the optimally adjusted random-count estimator. For a nonconstant population, the improvement over the unadjusted estimator is a factor independent of the budget: \[ \frac{\mathrm{Var}(\hat{x}^{\mathrm{random}})} {\mathrm{Var}(\hat{x}^{\mathrm{fixed}})} =1+\frac{\bar{x}^{\,2}}{\sigma^2}. \]
That factor can be large when the common offset is large compared with the spread of the values. In the baseline plot, its horizontal line touches the parabola at the unknown mean.
The correction we would have tuned is zero for every sample. Fixing the count has made the control variate implicit.
Without Replacement
If choosing distinct values is easy, or \(m\) is close to \(n\), it makes sense to sample without replacement. Draw a set \(S\) uniformly from the \(\binom nm\) subsets of exactly \(m\) distinct values and average them: \[ \hat{x}^{\mathrm{wor}}=\frac1m\sum_{i\in S}x_i. \]
Both without-replacement designs include each value with probability \(m/n\), giving unbiasedness [HT, JASA ’52]: \[ \mathbb{E}\!\left[\frac1m\sum_{i\in S}x_i\right] =\frac1m\sum_{i=1}^n x_i\frac mn =\bar{x}. \]
Claim: The fixed-count estimator without replacement has variance: \[ \mathrm{Var}(\hat{x}^{\mathrm{wor}}) =\frac{n-m}{n-1}\frac{\sigma^2}m. \]
Proof for the fixed-count sampler
Write \(\mathbf{1}[i\in S]\) for the indicator that value \(i\) is included. Because the size is fixed, the baseline cancellation lets us subtract \(\bar{x}\) inside the sum: \[ \hat{x}^{\mathrm{wor}}-\bar{x} =\frac1m\sum_{i=1}^n(x_i-\bar{x})\mathbf{1}[i\in S]. \]
We need the two indicator moments. Each indicator is Bernoulli with probability \(m/n\), so: \[ \mathrm{Var}(\mathbf{1}[i\in S]) =\frac mn\left(1-\frac mn\right). \]
For distinct \(i\) and \(j\), conditioning on \(i\) being included leaves \(m-1\) choices among \(n-1\) values: \[ \Pr(i,j\in S)=\frac mn\frac{m-1}{n-1}. \]
Subtracting the product of the marginal probabilities gives: \[\begin{align} \mathrm{Cov}(\mathbf{1}[i\in S],\mathbf{1}[j\in S]) &=\frac{m(m-1)}{n(n-1)}-\frac{m^2}{n^2}\\ &=\frac{m[(m-1)n-m(n-1)]}{n^2(n-1)}\\ &=-\frac{m(n-m)}{n^2(n-1)}. \end{align}\]
The second line puts the terms over a common denominator, and the third simplifies the bracket. The covariance is negative because including one value leaves fewer places for the others.
The centered values sum to zero. Squaring their sum and separating the diagonal terms therefore gives: \[\begin{align} \sum_{i\ne j}(x_i-\bar{x})(x_j-\bar{x}) &=\left(\sum_{i=1}^n(x_i-\bar{x})\right)^2-\sum_{i=1}^n(x_i-\bar{x})^2\\ &=-\sum_{i=1}^n(x_i-\bar{x})^2\\ &=-n\sigma^2. \end{align}\]
Expanding the variance of the weighted indicator sum now gives: \[\begin{align} \mathrm{Var}(\hat{x}^{\mathrm{wor}}) &=\frac1{m^2}\left[ \sum_i(x_i-\bar{x})^2\mathrm{Var}(\mathbf{1}[i\in S]) +\sum_{i\ne j}(x_i-\bar{x})(x_j-\bar{x}) \mathrm{Cov}(\mathbf{1}[i\in S],\mathbf{1}[j\in S]) \right]\\ &=\frac1{m^2}\left[ n\sigma^2\frac{m(n-m)}{n^2} +(-n\sigma^2)\left(-\frac{m(n-m)}{n^2(n-1)}\right) \right]\\ &=\frac{\sigma^2(n-m)}{mn}\left(1+\frac1{n-1}\right)\\ &=\frac{\sigma^2(n-m)}{mn}\frac n{n-1}\\ &=\frac{n-m}{n-1}\frac{\sigma^2}m. \end{align}\]
The second line substitutes the moments and the two sums of centered values. The third collects the shared factor, and the last two simplify it.
The finite population correction \((n-m)/(n-1)\) is at most one and reaches zero when we sample the whole population. Choosing distinct values gives an additional reduction beyond the one we already got by fixing the count.
For independent coin flips, let \(S\) contain each value independently with probability \(m/n\) and define: \[ \hat{x}^{\mathrm{coin}}=\frac1m\sum_{i\in S}x_i. \]
Claim: The random-count estimator without replacement has variance: \[ \mathrm{Var}(\hat{x}^{\mathrm{coin}}) =\left(1-\frac mn\right)\frac{\sigma^2+\bar{x}^{\,2}}m. \]
Proof for the random-count sampler
Here the indicators are independent, so their covariances are zero: \[\begin{align} \mathrm{Var}(\hat{x}^{\mathrm{coin}}) &=\frac1{m^2}\sum_{i=1}^n x_i^2\,\mathrm{Var}(\mathbf{1}[i\in S])\\ &=\frac1{m^2}\sum_{i=1}^n x_i^2\frac mn\left(1-\frac mn\right)\\ &=\frac{1-m/n}{m}\frac1n\sum_{i=1}^n x_i^2. \end{align}\]
The second line substitutes the Bernoulli variance and the third collects the constants. To recover the table’s expression, expand each raw square around the mean: \[\begin{align} \frac1n\sum_{i=1}^n x_i^2 &=\frac1n\sum_{i=1}^n\left[(x_i-\bar{x})^2+2\bar{x}(x_i-\bar{x})+\bar{x}^{\,2}\right]\\ &=\sigma^2+\frac{2\bar{x}}n\sum_{i=1}^n(x_i-\bar{x})+\bar{x}^{\,2}\\ &=\sigma^2+\bar{x}^{\,2}. \end{align}\]
The last line uses that the centered values sum to zero. Substituting gives: \[ \mathrm{Var}(\hat{x}^{\mathrm{coin}}) =\left(1-\frac mn\right)\frac{\sigma^2+\bar{x}^{\,2}}m. \]