2.1 Simulation and the Monte Carlo Method

Simulation will occupy most of our time in this module. The idea of drawing random variables from specified distributions and using the resulting ensemble of realisations to approximate quantities of interest is intoxicatingly powerful for a method which, at least in principle, is very simple.

Methods based around simulation are typically termed Monte Carlo² methods, following Metropolis and Ulam (1949). One characterization of such methods was given by Halton (1970):

Representing the solution of a problem as a parameter of a hypothetical population, and using a random sequence of numbers to construct a sample of the population, from which statistical estimates of the parameter can be obtained.

In some ways this is doing statistics backwards: rather than taking a real sample (of data) and attempting to infer parameters of an underlying population using analytical techniques, we devise a representation of a quantity of interest as a parameter of a hypothetical population and then obtain, artificially, a sample from that population before using the properties of this sample as a proxy for those of the population itself.

Perhaps an example might make it clear how simple this approach really is. You may well have seen this example before, but try to see it as a prototype for a general approach to approximate computation.

Example 2.1 (Computing \(\pi\) in the rain). Suppose we want to obtain an approximation of \(\pi\) using a simple experiment. Assume that we are able to produce “uniform rain” on the square extending to \(\pm1\) in two orthogonal directions, \([-1,1]^2 = [-1,1]\times[-1,1] = \{(x,y): x \in [-1,1], y \in [-1,1]\}\), such that the probability of a raindrop falling into a region \(\mathcal{R}\subseteq [-1,1]^2\) is proportional to the area of \(\mathcal{R}\), but independent of the position of \(\mathcal{R}\). It is easy to see that this is the case iff the two coordinates \(X,Y\) are independent realisations of uniform distributions on the interval \([-1,1]\) (in short \(X,Y \overset{\text{iid}}{\sim}{\textsf{U}{\left[-1,+1\right]}}\)).

Figure 2.1: Illustration of the estimation \(\pi\) using uniform rain.

Now consider the probability that a raindrop falls into the unit circle (see Figure 2.1). It is \[{\mathbb{P}}(\textrm{drop within circle})=\frac{\textrm{area of the unit circle}}{\textrm{area of the square}}=\frac{\underset{\{x^2+y^2\leq 1\}}{\int\!\int}1\;dxdy}{\underset{\{-1\leq x,y\leq 1\}}{\int\!\int}1\;dxdy}=\frac{\pi}{2\cdot 2}=\frac{\pi}{4}.\] In other words, \[\pi = 4 \cdot {\mathbb{P}}(\textrm{drop within circle}),\] i.e. there is an expression for the desired quantity \(\pi\) as a function of a probability.

Of course we cannot compute \({\mathbb{P}}(\textrm{drop within circle})\) without knowing \(\pi\). However, we can estimate the probability using our raindrop experiment. If we observe \(n\) raindrops, then the number of raindrops \(M\) that fall inside the circle is a binomial random variable, \(M \sim {\textsf{Bin}\left( n,p \right)}\) with \(p={\mathbb{P}}(\textrm{drop within circle})\).

Thus, if we observe that \(m\) raindrops fall within the circle, we approximate \(p\) using its maximum-likelihood estimate \(\hat p= {m}/{n},\) and we can estimate \(\pi\) by \(\hat \pi =4\hat p =4 \cdot \frac{m}{n}.\) Assume we have observed, as in Figure 2.1, that 77 of the 100 raindrops were inside the circle. In this case, our estimate of \(\pi\) is \(\hat \pi = {4\times 77}/{100}=3.08\), which is clearly some way from the truth.

However the strong law of large numbers (Theorem B.2) guarantees that the estimator \(\hat\pi = 4M/n\) converges almost surely to \(\pi\). Figure 2.2 shows the estimate obtained after \(n\) iterations as a function of \(n\) for \(n=1,\dots,2000\). You can see that the estimate improves as \(n\) increases.

Figure 2.2: Estimate of \(\pi\) together with approximate confidence intervals resulting from the raindrop experiment. Notice the jagged evolution of the confidence interval estimate: this is something which must be borne in mind when using simulation, if we use uncertainty estimates based upon estimated values then the quality of the estimate will determine also the quality of the uncertainty estimate. We can dramatically underestimate our own uncertainty if we are not careful.

We can assess the quality of our estimate by computing a confidence interval for \(\pi\). As we have \(X\sim {\textsf{Bin}\left( 100,p \right)}\), we can obtain a 95% confidence interval for \(p\) using a Normal approximation: \[\left[ 0.77 - 1.96 \cdot \sqrt{\frac{0.77\cdot (1-0.77)}{100}}, \; 0.77 + 1.96 \cdot \sqrt{\frac{0.77\cdot (1-0.77)}{100}}\right ] = [ 0.6875, \; 0.8525 ],\] As our estimate of \(\pi\) is four times the estimate of \(p\), we now also have a confidence interval for \(\pi\) which is simply \([2.750, \; 3.410]\).

In more general terms, let \(\hat \pi_n=4\hat p_n\) denote the estimate after having observed \(n\) raindrops. A \((1-2\alpha)\) confidence interval for \(\pi\) is \[\left[\hat \pi_n - z_{1-\alpha} \sqrt{\frac{\hat\pi_n(4-\hat\pi_n)}{n}}, \hat \pi_n + z_{1-\alpha} \sqrt{\frac{\hat\pi_n(4-\hat\pi_n)}{n}}\right].\]

Recall the main steps of this process:

• We have written the quantity of interest (in our case \(\pi\)) as an expectation (a probability is a special case of an expectation as \({\mathbb{P}}(A)={\mathbb{E}_{}\left[\mathbb{I}_A\right]}\) where \(\mathbb{I}_A(x)\) is the indicator function which takes value \(1\) if \(x\in A\) and 0 otherwise).

• We have replaced this algebraic representation of the quantity of interest with a sample approximation. The strong law of large numbers guaranteed that the sample approximation converges to the algebraic representation, and thus to the quantity of interest. Furthermore, the central limit theorem (Theorem B.3) allows us to assess the speed of convergence.

Of course we’d never use a method like this one to estimate \(\pi\); there are much faster ways of getting good estimates. Indeed, the rate of convergence here illustrates that these methods can be really computationally intensive. One major advantage of such methods, as we shall see, is that the rate of convergence of Monte Carlo estimates of expectations is independent of the dimension of the space on which the integral is defined, unlike more traditional approaches to numerical integration, and it is for hard problems in which other methods fail to produce any meaningful solution that simulation-based strategies are most useful.

2.2 Bootstrap Methods

Bootstrap methods are based around the, at first rather fanciful, idea that if we sample \(n\) times with replacement from a simple random sample of size \(n\) then the relationship between our resampled set of \(n\) points and the empirical distribution of the original sample is the same as the relationship between the original sample and its true distribution. This can be justified, at least asymptotically under regularity conditions via the Glivenko–Cantelli theorem (Theorem B.4) and its relatives. Sampling with replacement from a sample is equivalent to sampling from the distribution with its empirical distribution function.

There are many applications of bootstrap techniques, and many variations on the basic idea. In broad terms, the approach can be most useful when we want to construct confidence intervals without being able to compute the distribution of the test statistic. They can allow us:

• to mitigate certain types of bias; and

• to obtain similar results to those yielded by higher order asymptotic theory without doing the analysis that would be required to obtain them.

The paper of Efron (1979) is a reasonable first reference; a number of books on the subject have been written, including Efron (1982). If you already have Wasserman (2004) then you may find the informal introduction to these techniques presented in Chapter 8 helpful.

Bootstrap methods—and their strengths and limitations—will be discussed in the lectures.

2.3 Markov chains and Monte Carlo

Markov chains are objects with which you will become familiar during the APTS week, if you are not already. Roughly speaking, a Markov chain is a stochastic process for which the distribution of its future states is independent of its past given its current value. A more formal coverage of Markov chains was provided by the Applied Stochastic Processes module which took place in the previous APTS week.

The ergodic hypothesis was originally the work of Boltzmann and was intended to provide a characterisation of the long term average behaviour of a thermodynamic system. It said, approximately, that the time averaged behaviour of the microscopic configuration of a system was the same as the instantaneous average over a hypothetical ensemble of systems prepared in a particular way. In modern terms we would think of that hypothetical ensemble as a way of describing a probability distribution and we could think of the ergodic hypothesis as telling us (assuming it to be true) that the long-time-average behaviour of the stochastic process describing the evolution of the system coincided with an expectation with respect to that probability distribution.

The previous paragraph appears to be hinting at something like a law of large numbers and, indeed, that’s what a modern ergodic theorem describes. From a simulation point of view this is a tremendously powerful concept and there is an enormous literature on Markov chain Monte Carlo methods in which the trajectories of Markov chains are simulated for a long period of time and averages over these trajectories are calculated as a proxy for the expectation with respect to a particular distribution. A major contributing factor to the popularity and pre-eminence of these methods amongst computational statistics is that there exist recipes for the construction of Markov chains whose ergodic averages coincide with expectations with respect to essentially any distribution of interest, particularly the Metropolis (Metropolis et al. 1953) algorithm and its variants such as the Metropolis–Hastings algorithm (Hastings 1970) and the reversible jump algorithm (Green 1995). For some well-behaved statistical problems, the so-called Gibbs sampler (Geman and Geman 1984) provides what may seem an even more automatic approach, but it does come at the cost of some additional human calculation.

During the course of the APTS week we’ll see how to construct Markov chains to answer particular questions of interest and touch on all of the algorithms mentioned above.

There are vast numbers of resources which provide information on Markov chain Monte Carlo methods. I find that Robert and Casella (2004) covers this particular sub-field of Monte Carlo very well, but many other resources exist. For many years Gilks, Richardson, and Spiegelhalter (1996) has been something of a canonical reference, but has perhaps aged a little since its publication; perhaps less venerable but undoubtedly more up to date is Brooks et al. (2011).

2.4 Big Data is a Big Problem

Here’s the description of a dataset from a piece of work in which one of us was involved (Chan, Jenkins, and Song 2012):

We applied our method to samples from two populations of D. melanogaster (fruit flies): Raleigh, USA (RAL) and Gikongoro, Rwanda (RG). The RAL dataset consisted of the genomes of 37 inbred lines sequenced at a coverage of \(\geq 10\times\) by the Drosophila Population Genomics Project. The RG dataset comprised 22 genomes from haploid embryos sequenced at a coverage of \(\geq 25\times\) by the Drosophila Population Genomics Project 2…we were able to designate the ancestral allele in 1,755,040 of 2,475,674 high quality (quality score \(Q \geq 30\)) SNPs (single nucleotide polymorphisms) in the RAL sample (70.9%), and 2,213,312 out of 3,134,295 high quality SNPs in the RG sample (70.6%).

Each genome is summarised by its positions that differ from some other members of the sample (“SNPs”), so the data can be summarised by tables of size \(37 \times 2,475,674\) (RAL) and \(22 \times 3,134,295\) (RG). By the standards of modern statistics, these aren’t particularly large data sets, yet managing, storing and performing inference for data sets on this scale required considerable thought and effort.

The need for tools which can deal efficiently with data sets of this magnitude (and, indeed, much larger ones) is one of the reasons that computer intensive statistics is important. It would not be feasible to compute the sample mean of a data set with a quarter of a million elements without making use of a computer; of course, performing meaningful inference typically requires much more sophisticated computations than that.

A question to think about: from your personal perspective: how big is a large data set and how big is an enormous data set?

2.5 Warm-Up Exercises

A.1 Inference and Some Common Estimators

Computationally intensive methods are extremely prevalent in Bayesian statistics and people often make the mistake of thinking the two are synonymous. This is not the case. We will see in this module that computer intensive methods can be very useful in other areas of statistics, including likelihood-based inference. With this in mind, it is useful to recall a number of common estimation tasks we might want to carry out. It is very likely that you’ve come across all of these things before; the particular character of our interest is that we shall seek to find approximate solutions to these estimation problems in settings in which they are not analytically tractable (either because the computation is apparently not possible even in abstract terms or because carrying it out would take many times the age of the universe).

A.1.3 Hypothesis Tests

Closely related to the notion of a confidence interval is the hypothesis test. In order to distinguish between two possible explanations for observed data, a default scenario \(H_0\) termed the null hypothesis and an alternative \(H_1\), this procedure seeks to reject \(H_0\) when the data is in an appropriate sense unlikely to have arisen if \(H_0\) is true and not to reject \(H_0\) otherwise.

More precisely, given a test statistic \(T({\boldsymbol{X}})\), we seek a set of values \(C_\alpha\) such that \[\begin{aligned} {\mathbb{P}}(T({\boldsymbol{X}}) \in C_\alpha\mid H_0 \textrm{ is true}) &= \alpha, & {\mathbb{P}}(T({\boldsymbol{X}}) \in C_\alpha\mid H_1 \textrm{ is true}) > \alpha.\end{aligned}\] Often \(C_\alpha\) is obtained as the complement of an interval of values which are likely under the null hypothesis (viewed relatively, contrasting with the plausibility of those values under the particular alternative hypothesis). See Figure A.1 for an illustration.

Figure A.1: The critical or rejection region associated with a particular test; \(C_{\alpha} = \lbrack c,\infty)\) in this case.

Computing the distribution of \(T({\boldsymbol{X}})\) under \(H_0\) can be difficult (essentially impossible) in complicated situations. We will see that there are computational solutions to this problem.

A.2 Variability of Estimators and Uncertainty Quantification

A theme that turns out to be important in a number of forms in computer intensive statistics is the variability of estimate and the quantification of uncertainty. In classical statistics, the sampling distribution of an estimator is often used to provide some measure of uncertainty, perhaps via confidence intervals. The sampling distribution of a statistic \(T({\boldsymbol{X}})\) is simply the distribution which it has under repeated sampling of the data itself \({\boldsymbol{X}}\) from the model. In Bayesian statistics, the posterior distribution itself, \(p(\theta\mid {\boldsymbol{x}})\), summarises the uncertainty we have about the value of any parameter after incorporating the information contained in the data we have available. Both of these are things which we may wish to estimate using computer intensive statistics.

It is important to distinguish between the sampling distribution of a statistic or the posterior variance, which summarise in different ways the degree of uncertainty which must accompany any point estimate, and additional variability introduced by the procedure used to approximate an estimator. In particular, we will see that we often introduce additional auxiliary stochasticity during the computational procedures we use; this is undesirable and we seek to minimise it and to mitigate any influence it may have upon our estimation.

A.3 Warm-Up Exercises

The main purpose of the following is to show some places in which it quickly becomes difficult to proceed analytically, but for which computational methods might work well, and to highlight the types of quantities which we will want to be able to compute or approximate in this module.

B.1 A Word on Probability

We won’t see technical (measure theoretic) probability in this module—we can manage without it for our purposes. There are a few places where this may cause us a slightly loss of generality, but these will be highlighted.

When we talk about probability we assume that there is some underlying sample space, \(\Omega\), from which exactly one outcome occurs every time our experiment is realised: an elementary event, \(\omega \in \Omega\). A random variable can be thought of as a measurement which we could make when the experiment is carried out and can be modelled mathematically as a function which maps the sample space to the real numbers, \(X : \Omega \rightarrow \mathbb{R}\) (see Figure B.1).

Figure B.1: Random variables as functions from the sample space, \(\Omega\), to \(\mathbb{R}\). If \(\omega\) occurs in an experiment then \(X\) will be measured as having the value \(X(\omega)\). If we want to know whether \(X\) takes a value in the set \(A\) (the green rectangle) then we need to check whether \(\omega\) takes a value which \(X\) maps into \(A\), this set is often denoted \(X^{-1}(A)\) and corresponds to the blue ellipse here.

When we talk about the probability of some event \(B \subseteq \Omega\) we mean exactly \({\mathbb{P}}(B) := {\mathbb{P}}(\omega \in B)\), i.e. the probability that the elementary outcome which occurs is contained within \(B\). When we talk about the probability of a random variable, \(X\), taking a value in a set, e.g. \({\mathbb{P}}(X \in A)\) we really mean the probability that the elementary outcome that occurs is such that \(X\) takes a value in \(A\): \({\mathbb{P}}(X \in A) := {\mathbb{P}}(\{\omega:X(\omega) \in A\}) = {\mathbb{P}}(\omega \in X^{-1}(A))\) where \(X^{-1}(A)\) denotes the pre-image of \(A\) under \(X\), i.e. the collection of all points in \(\Omega\) which are mapped into \(A\) by \(X\): \(X^{-1}(A) = \{\omega:X(w) \in A\}\).

B.2 Convergence In Probability

A sequence of random variables \(X_1,X_2,\dots\) is said to converge in probability to a limiting value \(x\), written \(\lim_{n\rightarrow\infty} X_n = x \text{ (in probability)}\) or \(X_n \overset{p}{\rightarrow}x\), if for every \(\epsilon > 0\) the probability that \(X_n\) is further from \(x\) than \(\epsilon\) converges to zero, i.e. \[\lim_{n \rightarrow \infty} {\mathbb{P}}\left(|X_n - x| > \epsilon \right) = 0.\]

The first theoretical result of interest tells us that in this sense, averages obtained from simple samples will converge to the population average.

B.3 Almost Sure Convergence

A sequence of random variables \(X_1,X_2,\dots\) is said to converge almost surely (a.s.) or converge with probability 1 to a limiting value \(x\), written \(\lim_{n\rightarrow\infty} X_n = x \text{ (a.s.)}\) or \(X_n \overset{a.s.}{\rightarrow}x\), if \({\mathbb{P}}\left(X_n \rightarrow x \right)={\mathbb{P}}\left(\left\{ \omega : X_n(\omega) \rightarrow x \right\}\right) = 1\).

Almost sure convergence is strictly stronger than convergence in probability. We can, however, be sure that under weak assumptions the average from a simple random sample will convergence almost surely to the underlying population average.

B.4 Some Ideas Related to Convergence of Distributions

B.5 Warm-Up Exercises

If you feel like you could do with reminding yourself how these ideas of stochastic convergence work then you might like to have a go at the following; if you already know how to answer them then don’t waste your time and, similarly, if you have no interest in such things then it won’t be essential to the lectures for this module that you have worked these out.