Random tricks for computing costly sums

The uniform estimators discussed above are a good start, but they can fall short of being precise. This happens in particular in "needle in a haystack" situations where the sum is highly imbalanced, the vast majority of the terms \(f(x)\) contributing negligibly to the sum \(S\). Then the variance is high, and a very high sample size is needed to achieve good precision.

In such cases, it pays to oversample the terms \(x\) that are most likely to contribute significantly to the sum. This is the principle of Importance Sampling, in which the \(x\) are drawn with different selection probabilities, and the sum is weighted to compensate for the bias this introduces.

Concretely, we draw the sample \(\left(X_j\right)_j\) i.i.d. but not uniformly: each \(x \in \mathcal{X}\) has a probability \(p(x)\) of being selected (\(p(\cdot)\) is called the proposal distribution). The estimator is then:

\begin{equation} \label{org5661861} \hat{S}_p := \frac{1}{n} \sum_{j}{\frac{1}{p(X_j)} f(X_j)} \end{equation}

This is also an unbiased estimator: \(\hat{S}_p \overset{\mathbb{E}}{\approx} S\), which is easily verified as follows:

\begin{align*} \mathbb{E}\left[\hat{S}_p\right] = \frac{n}{n} \mathbb{E}\left[\frac{1}{p(X_1)} f(X_1)\right] = \sum_{x}{p(x) \frac{1}{p(x)} f(x)} = S \end{align*}

What can make \(\hat{S}_p\) more interesting than \(\hat{S}\) is its precision characteristics. The lowest possible variance is achieved when \(p(x)\) is proportional to \(|f(x)|\). This optimum is typically impossible to achieve (you'd have to compute \(\sum_x{|f(x)|}\)) but usually you can rely on some heuristic pointing you to a small subset of \(\mathcal{X}\) where most of the sum is concentrated.

As before, other sampling schemes can be used to draw the sample \(\left(X_j\right)_j\). The proposal distribution \(p(x)\) can be used to draw without replacement. Bernoulli sampling can also be used with some non-constant selection probability function \(q(x) \in (0; 1]\), and the formula becomes:

\begin{equation} \label{orgffa4bf4} \hat{S}_B := \sum_{j}{\frac{1}{q(X_j)} f(X_j)} \end{equation}

Equation \eqref{orgffa4bf4} also works in an Poisson sampling scheme, in which the number of repetitions of \(x\) in the sample follows \(\text{Poisson}(\lambda = q(x))\), in which case \(q(x) > 1\) is allowed.

Anticipating the next sections, we encourage readers to think of \(q(x)\) and \(n p(x)\) as the expected number of times that \(x\) is selected.

Table 2 describe the variance of these estimators, and how it can be estimated from the sample:

Table 2: Variance of estimators from various weighted sampling schemes. The rightmost column is an unbiased estimator of \(\text{Var}\left[\hat{S}\right]\).

Sampling scheme	\(\text{Var}\left[\hat{S}\right]\)	Variance estimator
With replacement	\(\frac{1}{n} \text{Var}\left[ \frac{1}{p(X_1)}f(X_1) \right]\)	\(\frac{1}{n-1} \left( \left(\sum_j{\frac{1}{p(X_j)}f(X_j)^2} \right) - \hat{S}^2 \right)\)
Without replacement	\(\frac{1}{n} \text{Var}\left[ \frac{1}{p(X_1)}f(X_1) \right] \left(1 - \frac{n-1}{\lvert X \rvert - 1}\right)\)	\(\frac{1}{n-1} \left( \left(\sum_j{\frac{1}{p(X_j)}f(X_j)^2} \right) - \hat{S}^2 \right) \left(1 - \frac{n-1}{\lvert X \rvert - 1}\right)\)
Bernoulli	\(\sum_x{\frac{1 - q(x)}{q(x)} f(x)^2}\)	\(\sum_j{\frac{1 - q(X_j)}{q(X_j)^2} f(X_j)^2}\)
Poisson	\(\sum_x{\frac{1}{q(x)} f(x)^2}\)	\(\sum_j{\frac{1}{q(X_j)^2} f(X_j)^2}\)

Stratified estimation is an estimation technique with similar characteristics to Importance Sampling: we sample more densely in some parts of the space than others. The idea is to partition \(\mathcal{X}\) into \(K\) disjoint strata \(\mathcal{X}_1, \mathcal{X}_2, \dots, \mathcal{X}_K\):

\begin{equation} \label{org71824a1} \mathcal{X} = \mathcal{X}_1 \uplus \mathcal{X}_2 \uplus \dots \uplus \mathcal{X}_K \end{equation}

In each stratum \(\mathcal{X}_k\), we draw a sample \(\left(X_{kj}\right)_{j=1}^{n_k}\) of size \(n_k\). Typically, we'll have at least one stratum \(k\) which we expect to contain most of the large \(f(x)\) terms, and sample abundantly from it even though it makes from a small subset of \(\mathcal{X}\).

The corresponding estimator is then:

\begin{equation} \label{org95735c6} \hat{S}_\text{strat} := \sum_{k=1}^{K}{\frac{|\mathcal{X}_k|}{n_k} \sum_{j=1}^{n_k}{f(X_{kj})}} \end{equation}

Because strata are sampled independently, the variance of \(\hat{S}_\text{strat}\) is simply the sum of the per-stratum variances:

\begin{equation} \label{org7266da9} \text{Var}\left[\hat{S}_\text{strat}\right] = \sum_{k=1}^{K}{\text{Var}\left[\hat{S}_k\right]} \end{equation}

Stratified estimation can be further refined by using the ideas of previous sections within each stratum. At which point, our list of estimators is starting to look like a zoo. It's time for some unification.

We have seen how to estimate \(S\) using a small sub-sample of terms, and how to make the most of the sample size using weighted sampling schemes. How could we possibly want more? Well, here's one situation where we could want more. Suppose that you used Importance Sampling with sample size \(N_1\) and proposal distribution \(p_1(x)\). Then your colleague Meryem tells you that she used Bernoulli sampling with selection probability \(q_2(x)\). It would be a shame not to combine your samples to get better precision, wouldn't it? Spoiler alert, here's how to do it:

\begin{equation} \hat{S}_{1 \cup 2} := \sum_j{\frac{1}{N_1 p_1(X_{j}) + q_2(X_{j})}f(X_{j})} \end{equation}

in which the sum iterates over your merged samples. I will not bother to prove that this works in this particular use case - we will soon see a general result that will make this sort of formula obvious.

As another example, imagine that you are dealing with high-resolution geospatial data spanning all of India. Your data consists of rasters with 30m pixels: \(\mathcal{X}\) is the set of pixels and \(f(x)\) consists of computing the average distance of \(x\) to an urban area. No single file can contain so many pixels, and so the data is split in 60km square tiles, one raster file per tile.

This is an example where the cost structure of your computation is too complicated for the previous estimators to be useful. You cannot just scatter the sample all over the landscape. Instead, you'll want to use a multi-stage sampling scheme: first sample a subset of tiles, and then maybe sub-sample the pixels within each tile.

This time I won't show you the formula - can you come up with an estimator that works?

The first step towards unifying the previous results is to notice that all the estimators we studied above can be written in the form:

\begin{equation} \label{org0e22f27} S \overset{\mathbb{E}}{\approx} \hat{S} = \sum_{j}{\frac{1}{\nu(X_j)} f(X_j)} \end{equation}

in which \(\nu(x)\) is the expected number of times that \(x\) is selected in the sample:

\begin{equation} \label{orgfbd389f} \nu(x) := \mathbb{E}\left[\hat{n}_x\right] \end{equation}

\(\nu(\cdot)\) might be described as the "expected count" function; I will call it the expected frequency mass function.

The proof to equation \eqref{org0e22f27} is elementary:

\begin{align*} \mathbb{E}\left[\hat{S}\right] &= \mathbb{E}\left[\sum_{x}{\sum_{j | X_j = x}{\frac{1}{\nu(X_j)} f(X_j)}}\right] \\ &= \mathbb{E}\left[\sum_{x}{\frac{1}{\nu(x)} f(x) \sum_{j | X_j = x}{1}}\right] \\ &= \mathbb{E}\left[\sum_{x}{\frac{1}{\nu(x)} f(x) \hat{n}_x}\right] \\ &= \sum_{x}{\frac{1}{\nu(x)} f(x) \mathbb{E}\left[\hat{n}_x\right]} \\ &= \sum_{x}{\frac{1}{\nu(x)} f(x) \nu(x)} \\ &= \sum_{x}{f(x)} \\ \mathbb{E}\left[\hat{S}\right] &= S \end{align*}

Equation \eqref{org0e22f27} gives us a general and versatile method for estimating sums:

1. Design a sub-sampling scheme that fits your problem. Your intuition probably has good approaches to suggest, and you can also use some of the techniques listed above, like Importance Sampling or Stratified Estimation.
2. Find the expected frequency mass function \(\nu(x)\) that corresponds to your sub-sampling scheme.
3. Draw a random sample \(\left(X_j\right)_j\) from your sub-sampling scheme.
4. Compute and record \(\left(f(X_j)\right)_j\).
5. Estimate the sum by applying equation \eqref{org0e22f27}.

Unfortunately, there is no general formula for the variance of the estimator defined in equation \eqref{org0e22f27}. The variance will depend on the specifics of the sampling procedure; it's possible to contrive sampling procedures of arbitrarily bad variance, by not sampling densely in the important parts of the space and/or by making samples highly correlated. The variance will have to be worked out specifically for each sampling procedure, which can usually be done by applying basic mathematical properties of variance: variance responds quadratically to scaling (\(\text{Var}[\lambda Y] = \lambda^2 \text{Var}[Y]\)), is additive on independent random variables (\(\text{Var}[Y_1 + Y_2] = \text{Var}[Y_1] + \text{Var}[Y_2]\)), can be expressed from moments (\(\text{Var}[Y] = \mathbb{E}[Y^2] - \mathbb{E}[Y]^2 = \mathbb{E}[(Y - \mathbb{E}[Y])^2]\)).

In addition, the conditional variance formula is useful for multi-stage sampling schemes:

\begin{equation} \label{org093c7b6} \text{Var}[Y] = \mathbb{E}\left[\text{Var}[Y | Z]\right] + \text{Var}\left[\mathbb{E}[Y | Z]\right] \end{equation}

Finally, the variance of a sum with correlated terms can be worked out from the covariances of the terms:

\begin{equation} \label{org4943892} \text{Var}\left[ \sum_{i=1}^m{Y_i} \right] = \sum_{i,j}{\text{Cov}\left[ Y_i, Y_j \right]} \end{equation}

We will now make the above results applicable to a wider range of problems, like estimating continuous integrals and probabilistic expectations.

We will need some notions from Measure Theory; this is unsurpring, since Measure Theory is essentially the mathematical theory of sums. Said informally, a measure is something that distributes of positive mass over a space \(\mathcal{X}\). Concretely, a measure \(\mu\) is a function that assigns a non-negative "mass" \(\mu(A)\) to each set \(A \subseteq \mathcal{X}\) (with some constraints of what \(A\) and \(\mu(A)\) are allowed to be, such as non-negativity, addivity, etc.). Furthermore, by definition of integrals, we have:

\begin{equation} \label{orgadf7ccf} \int_{x \in A}{d\mu(x)} = \mu(A) \end{equation}

For example, if \(\mathcal{X} = \left\{x_1, x_2, \dots \right\}\) is a discrete set, the counting measure \(\kappa\) over \(\mathcal{X}\) maps each set \(A \subseteq \mathcal{X}\) to the number of elements \(x_i \in A\):

\begin{equation} \label{org6d0d23b} \kappa(A) := |A| \end{equation}

As such, our discrete sum \(S\) (defined by equation \eqref{org2fba7dc}) can be rewritten as an abstract integral:

\begin{equation} \label{org7317578} S := \sum_{x \in \mathcal{X}}{f(x)} = \int_{x \in \mathcal{X}}{f(x) d\kappa(x)} \end{equation}

This more general framing will prove useful when we define the estimand \(S\) to be something else than a discrete sum.

Another example of measure is the Dirac distribution \(\delta_x\) at a point \(x \in \mathcal{X}\), which maps each set \(A \subseteq \mathcal{X}\) to \(1\) if \(x \in A\) and \(0\) otherwise:

\begin{equation} \label{orga49ab20} \delta_{x}(A) := \begin{cases} 1 & \mbox{if } x \in A \\ 0 & \mbox{if } x \notin A \end{cases} \end{equation}

Measure theory tells us that an integral over a measure \(\mu\) can be rewritten as an integral over another measure \(\nu\), via the integral reweighting trick:

\begin{equation} \label{org5db7389} \int_{x \in \mathcal{X}}{g(x) d\mu(x)} = \int_{x \in \mathcal{X}}{g(x) \frac{d\mu}{d\nu}(x) d\nu(x)} \end{equation}

in which \(\frac{d\mu}{d\nu}(x)\) is the Radon-Nikodym derivative of \(\mu\) with respect to \(\nu\). No, don't run, there's an intuitive concept behind this scary name. Said informally, \(\frac{d\mu}{d\nu}(x)\) is a concentration ratio: it says by what factor \(x\) is more weighted by \(\mu\) than by \(\nu\). Equation \eqref{org5db7389} only works if \(\nu\) dominates \(\mu\), which essentially means that \(\frac{d\mu}{d\nu}\) is well defined: formally, it means that \(\mu(A) = 0\) whenever \(\nu(A) = 0\).

We now apply this theory to the use of random samples. So far, we had seen our random sample \((X_j)_j\) as a point process over \(\mathcal{X}\), that is, a random set of points in \(\mathcal{X}\) (more rigorously, a random collection of points, since duplicates are allowed). We now see our random sample as a random measure \(\hat{n}_X\) over \(\mathcal{X}\), which maps each set \(A \subseteq \mathcal{X}\) to the number of sample points \(X_j\) "landing" inside \(A\):

\begin{equation} \label{org2a2a3fd} \hat{n}_X(A) := \sum_{j | X_j \in A}{1} = \sum_{j}{\delta_{X_j}(A)} \end{equation}

Again, this measure-theoretic framing lets us rewrite discrete sums over the sample as abstract integrals over measure \(\hat{n}_X\):

\begin{equation} \label{org83b36a9} \sum_j{g(X_j)} = \int_{x \in \mathcal{X}}{g(x) d\hat{n}_X(x)} \end{equation}

The notion of probabilistic expectation can be applied to random measures. The expectation measure of \(\hat{n}_X\) is the measure \(\nu\) defined as:

\begin{equation} \label{org2da55f9} \nu(A) := \mathbb{E}\left[\hat{n}_X(A)\right] \end{equation}

Then, probability theory allows us to exchange the integral and expectation operators:

\begin{equation} \label{orgdfab63d} \mathbb{E}\left[\int_{x \in \mathcal{X}}{g(x) d\hat{n}_X(x)}\right] = \int_{x \in \mathcal{X}}{g(x) d\nu(x)} \end{equation}

By the way, this "point process viewpoint" can make the Radon-Nikodym derivative \(\frac{d\mu}{d\nu}\) more intuitive. If \((X_{j})_j\) and \((X'_{j})_j\) are two point processes with expectation measures \(\nu\) and \(\nu'\), then the Radon-Nikodym derivative \(\frac{d\nu'}{d\nu}(x)\) tells us how much more often \(x\) gets drawn in \(X'\) than in \(X\); or said differently, it tells us by what factor \(X'\) over-represents \(x\) compared to \(X\).

We now have all the tools we need to estimate general integrals using random samples. We no longer assume that \(\mathcal{X}\) is a discrete space. Instead, we assume that our estimand \(S\) is now an integral over a measure \(\mu\):

\begin{equation} \label{org12799a5} S := \int_{x \in \mathcal{X}}{f(x) d\mu(x)} \end{equation}

We now draw a sample \((X_j)_j\) with an expectation measure \(\nu\) that dominates \(\mu\). By combining equations \eqref{org83b36a9}, \eqref{orgdfab63d} and \eqref{org5db7389}, we obtain an unbiased estimator of \(S\):

\begin{align*} S &= \int_{x \in \mathcal{X}}{f(x) d\mu(x)} \\ &= \int_{x \in \mathcal{X}}{f(x) \frac{d\mu}{d\nu}(x) d\nu(x)} \\ &= \mathbb{E}\left[ \int_{x \in \mathcal{X}}{\frac{d\mu}{d\nu}(x) f(x) d\hat{n}_X(x)} \right] \\ S &= \mathbb{E}\left[ \sum_j{\frac{d\mu}{d\nu}(X_j) f(X_j)} \right] \end{align*}

In summary, we estimate \(S\) through the sample-reweighting estimator \(\hat{S}_X\):

\begin{equation} \label{org77807e8} S \overset{\mathbb{E}}{\approx} \hat{S}_X := \sum_j{\frac{d\mu}{d\nu}(X_j) f(X_j)} \end{equation}

Taking a step back, we can say that we are essentially approximating the measure \(\mu\) using a weighted mixture of Dirac distributions:

\begin{equation} \label{org594a7de} \mu \overset{\mathbb{E}}{\approx} \sum_j{\frac{d\mu}{d\nu}(X_j) \delta_{X_j}} \end{equation}

Let us now revisit the discrete case (equation \eqref{org0e22f27}) from this more general perspective. As we saw in equation \eqref{org7317578}, \(\mu\) is simply the counting measure \(\kappa\) over \(\mathcal{X}\), with mass function \(\kappa(x) = 1\); the expectation measure \(\nu\) has mass function \(\nu(x)\), so the Radon-Nikodym derivative \(\frac{d\mu}{d\nu}\) is concretely expressed as:

\begin{equation} \label{orgd429edf} \frac{d\mu}{d\nu}(x) = \frac{1}{\nu(x)} \end{equation}

from which equation \eqref{org0e22f27} follows directly as a special case of equation \eqref{org77807e8}.

A common use case is to estimate some expected value. We now assume that the estimand \(S\) is the expected value of \(f(X)\) in which \(X\) follows a probability distribution with density \(p(x)\):

\begin{equation} S := \mathbb{E}\left[f(X)\right] = \int_{x \in \mathcal{X}}{f(x) p(x) dx} \end{equation}

The most common to estimate \(S\) is to draw an i.i.d. sample \((X_j)_{j=1}^n\) from \(p\), and compute the sample mean:

\begin{equation} \label{org4c9e3f1} S \overset{\mathbb{E}}{\approx} \hat{S}_\frac{1}{n} := \frac{1}{n} \sum_{j=1}^n{f(X_j)} \end{equation}

It's interesting to revisit the sample mean as a special case of our more general framework. Our i.i.d. sampling procedure can be viewed as a point process; its expectation measure is \(\nu = n p\), with density function \(n p(x)\). The Radon-Nikodym derivative is:

\begin{align*} \frac{dp}{d\nu}(x) = \frac{dp/dx}{n dp/dx} = \frac{p(x)}{n p(x)} = \frac{1}{n} \end{align*}

so that the sample mean estimator (equation \eqref{org4c9e3f1}) follows directly as a special case of the sample-reweighting estimator (equation \eqref{org77807e8}):

\begin{equation} \hat{S}_X := \sum_j{\frac{1}{n} f(X_j)} = \frac{1}{n} \sum_{j=1}^n{f(X_j)} = \hat{S}_\frac{1}{n} \end{equation}

Some sampling procedures like Markov Chain Monte Carlo (MCMC) do not draw independent samples (the chains have autocorrelation). However, once ergodicity is achieved, a chain of length \(n\) does have \(n p\) as its expectation measure, so that equation \eqref{org4c9e3f1} still applies (with somewhat degraded variance).

Now assume that the sample \((X_j)_{j=1}^n\) was drawn from some other probability distribution \(p'\), with density function \(p'(x)\). Then applying equation \eqref{org77807e8} gives us what has been called in the literature the likelihood ratio estimator:

\begin{equation} \label{org4d0939e} S \overset{\mathbb{E}}{\approx} \hat{S}_\frac{p}{n p'} := \frac{1}{n} \sum_{j=1}^n{\frac{p(x)}{p'(x)} f(X_j)} \end{equation}

The above is how Importance Sampling is generally framed in the literature. Equation \eqref{org4d0939e} also works in the discrete case, i.e. when \(p(x)\) and \(p'(x)\) are probability mass functions instead of probability density functions.

Recall that the proof of equation \eqref{org77807e8} involved the expectation of an integral over the sample viewed as a random measure \(\hat{n}_X\):

\begin{equation} \label{org3cb0833} S \overset{\mathbb{E}}{\approx} \int_{x \in \mathcal{X}}{\frac{d\mu}{d\nu}(x) f(x) d\hat{n}_X(x)} \end{equation}

We will now simply forget that \(\hat{n}_X\) represents a point sample, thus making equation \eqref{org3cb0833} more broadly applicable. This is interesting because sometimes, we have more efficient ways of computing integrals than by drawing a point sample.

In the United States, most of the damage caused by wildfires to structures is due to firebrands - burning embers that get carried by the wind and ignite buildings. The risk caused by firebrands can be estimated by running Monte Carlo simulations that model hypothetical fires and their firebrands emissions.

The most basic way in which this is done it to simulate a large number of fires \((F_j)_j\), and for each fire to simulate the trajectories and structure ignitions of firebrands emitted by burning fuels, recording which structures end up ignited. Effectively, we're estimating structure ignition probability by drawing a point process \((B_j(s))_j\) of structure ignitions, \(B_j(s)\) being a Bernoulli random variable corresponding to \(j\) burning \(s\), and the probability² \(\beta(s)\) of structure \(s\) burning across all potential fires is estimated as:

\begin{equation} \beta(s) = \int_{F}{\mathbb{P}\left[s \text{ burns} | F\right] d\mu(F)} \overset{\mathbb{E}}{\approx} \hat{\beta}(s) := \sum_j{\frac{d\mu}{d\nu}(F_j) B_j(s)} \end{equation}

in which \(\mu\) and \(\nu\) are respectively the expectation measures of real-world and simulated fire occurrence (which are known by design of the simulations).

This estimation could be done more efficiently, however. Oftentimes, the stochastic model for firebrand trajectories is simple enough to be analytically tractable: typically, each burning pixel emits a random number of firebrands based on its burning conditions, then each firebrand "jumps" independently, with a log-normal jump in the wind direction and a normal jump orthogonal to it, so that we can compute the spatial probability density of the jump. Then we might apply a probability of successful structure ignition based on things like travelled distance and structure characteristics. The upshot is: for each structure \(s\) and each burned pixel \(z\), we know how to compute the expected number \(\phi_j(s;z)\) of firebrands originating from \(z\) that burn \(s\) conditional on fire \(j\).

Then the total expected number \(\phi_j(s)\) of firebrands igniting structure \(s\) during fire \(j\) is simply the sum over burned pixels:

\begin{equation} \phi_j(s) = \sum_{z}{\phi_j(s;z)} \end{equation}

It's sensible to assume that the number of firebrands actually burning \(s\) is Poisson-distributed, and therefore the probability that the structure burns during fire \(j\) is \(\beta_j(s) = \left(1 - e^{- \phi_j(s)}\right)\). Then the aggregated risk estimate becomes:

\begin{equation} \beta(s) \overset{\mathbb{E}}{\approx} \sum_j{\frac{d\mu}{d\nu}(F_j) \beta_j(s)} \end{equation}

Of course, all of the techniques that we have seen in previous sections are applicable to sampling the simulated fires. For example, we might use Importance Sampling to oversample the fires that are most likely to emit firebrands, for example because they happen in firebrand-prone vegetation like eucalyptus.