9.1 Introduction

Computer simulation, also known as stochastic simulation or Monte Carlo simulation, is a virtual experiment in which we mimic the biological process on a computer to study its properties. It is particularly useful for studying complex systems or problems which are analytically intractable. Use of random numbers is a major feature of the method. Some authors use the term Monte Carlo simulation when the answer is deterministic, as in calculation of an integral by Monte Carlo integration, and stochastic simulation or computer simulation when the answer involves random variations, as in studies of biases and variances of an estimation method. We ignore such a distinction here.

Simulation is a useful approach for validating a theory or program implementation when the method of analysis is complex. When the model is intractable analytically, simulation provides a powerful way of studying it. Indeed, as it is easy to simulate even under a complex model, there is no need to make overly simplistic assumptions. Simulation is commonly used to compare different analytical methods, especially their robustness when the underlying assumptions are violated, in which case theoretical results are often unavailable. Simulation also plays a useful role in education. By simulating under a model and observing its behaviour, one gains intuitions about the system. Lastly, Monte Carlo simulation forms the basis of many modern computation-intensive statistical methods, such as bootstrapping, importance sampling, and Markov-chain Monte Carlo (see Chapter 5).

A consequence of the ease of simulating and the availability of cheap and fast computers is that simulation is widely used and misused. Many simulation studies are poorly designed and analysed, with results largely uninterpretable. One should bear in mind that simulation is experimentation, and the same care should be exercised in the design and analysis of a simulation experiment as an ordinary experiment. While it is perhaps the simplest thing one can do with a model, simulation does offer the uninitiated ample opportunities for mistakes. A major limitation of simulation is that one can examine only a small portion of the parameter space, and the behaviour of the system may be different in other unexplored regions of the parameter space. One should thus resist the temptation of over-generalization. Analytical results are in general superior as they apply to all parameter values.

This chapter provides a cursory introduction to simulation techniques. The reader is invited to implement ideas discussed here. Any programming language, such as (p. 294 ) C/C++, Fortran, BASIC, Java, Perl, or Python, will serve the purpose well. The reader may refer to many of the textbooks on simulation, such as Ripley (1987) and Ross (1997).

9.2 Random number generator

Random variables from the uniform distribution U (0, 1), called random numbers, are fundamentally important in computer simulation. For example, to simulate an event (E) that occurs with probability 0.23, we draw a random number u ∼ U (0, 1). If u < 0.23, we decide that the event occurs; otherwise it does not (Fig. 9.1). It is obvious that the probability that u < 0.23 is 0.23. The U (0, 1) random numbers form the basis for generating random variables from other distributions.

A mathematical algorithm is used to generate a sequence of numbers between 0 and 1 that appear random. The algorithm is deterministic and will produce the same fixed sequence of numbers every time it is run. The generated numbers are thus not random at all, and are more appropriately called pseudo-random numbers but often simply random numbers. The mathematical algorithm is called a random-number generator. A typical simulation study may need millions or billions of random numbers, so it may be important to have a reliable and efficient random-number generator.

A commonly used algorithm is the multiplication-congruent method, given by the following formulae (9.1) (9.2) Here, A_i, c, and M are all positive integers: c is called the multiplier and M the modulus, with (a mod M) to be the remainder when a is divided by M. cA_i−1 and A_i have the same remainder when divided by M and are said to be congruent. A₀ is the initial value, known as the seed. A₁, A₂, … will be a sequence of integers between 0 and M – 1, and u₁, u₂, … will be a sequence of (pseudo-) random numbers.

Fig. 9.1 Sampling from a discrete distribution. (a) To simulate an event that occurs with probability 0.23, draw a random number u ∼ U (0, 1). If it falls into the interval (0, 0.23), the event occurs (E). Otherwise the event does not occur (Ē). (b) Similarly to sample a nucleotide with probabilities 0.1, 0.2, 0.3, and 0.4, for T, C, A, and G, respectively, draw a random number u ∼ U (0, 1) and choose the corresponding nucleotide depending on whether u falls into the four intervals of lengths 0.1, 0.2, 0.3, and 0.4.

(p. 295 ) As an example, imagine a computer using the familiar decimal system, with the register (storage) having only four digits. The computer can thus represent any integer in the range (0000, 9999). Let M = 10⁴, c = 13, and A₀ = 123. Equation (9.1) then generates the A_i sequence 1599, 6677, 1271, 6333, 8959, …, corresponding to the u_i sequence 0.1599, 0.6677, 0.1271, 0.6333, 0.8959, …. Note that if the product cA_i−1 is greater than 9999, we simply ignore the high digits to get the remainder A_i; most computer systems allow such ‘overflow’ without any warning. The sequence u_i does not show an obvious trend and appears random. However, as this computer can represent only 10 000 integers and many of them will never occur in the sequence, very soon a certain number A_i will reoccur and then the sequence will repeat itself. A real computer, however, has more digits, and by intelligent choices of M and c (and, to a lesser extent, A₀ as well), the period of the sequence can be made very long.

Real computers use the binary system. Thus the natural choice of M is 2^d, often with d to be the number of bits in an integer (say, 31, 32, or 64). Note that such a choice of M eliminates the need for the division in equation (9.1) for the modulus calculation and also makes the division in equation (9.2) a trivial calculation.

How do we decide whether the random numbers generated are random enough. The requirement is that the numbers generated by the algorithm should be indistinguishable from random draws from the uniform distribution U (0, 1). This can be examined using statistical tests, examining different departures from the expectation. For example, the random numbers should have the correct mean (1/2) and variance (1/12), and should not be autocorrelated. In general, it is not advisable for us biologists to devise random-number generators. Instead one should use algorithms that have been well tested (Knuth 1997; Ripley 1987). Ripley (1987) discussed various tests and commented that random-number generators supplied by computer manufacturers or compiler writers were not trustable. It is unclear whether the situation has improved in the last two decades. For the exercises of this chapter, it is sufficient to use a generator provided by the programming language.

9.3 Generation of continuous random variables

Random variables used in a computer simulation are all generated using U (0, 1) random numbers. A very useful method is the transformation method, based on the fact that a function of a random variable is itself a random variable (see Appendix A). Thus we generate a random number u ∼ U (0, 1) and apply an appropriate variable transform, x = g(u), to generate a random variable x with the desired distribution.

An important transformation method is the inversion method. Suppose random variable x has a cumulative density function (CDF) F(x). Consider u = F(x) as a function of x so that u is itself a random variable. Its distribution is known to be uniform: u = F(x) ∼ U (0, 1). Thus, if the inverse transform x = F⁻¹(u) is available analytically, one can use it to generate x using uniform random number u.

(p. 296 )

1. 1. Uniform distribution. To generate a random variable from the uniform distribution x ∼ U (a, b), one generates u ∼ U (0, 1) and apply the transform x = a + u(b − a).

2. 2. Exponential distribution. The exponential distribution with mean θ has density function f (x) = θ⁻¹e^−x/θ and CDF $F(x)={\int }_0^x{\text{e}}^{-y/\theta }\text{d}\text{y}=1-{\text{e}}^{-x/\theta }$. Let u = F(x). This function is easy to invert, to give x = F⁻¹(u) = − θ log (1 − u). Thus we generate u ∼ U (0, 1) and then x = − θ log (1 − u) will have the exponential distribution with mean θ. Since 1 − u is also U (0, 1), it suffices to use x = −θ log (u).

3. 3. Normal distribution. Very often the CDF is not easy to invert. This is the case with the CDF of the normal distribution, which can only be calculated numerically. However, Box and Muller (1958) described a clever method for generating standard normal variables using variable transform. Let u₁ and u₂ be two independent random numbers from U (0, 1). then (9.3) are two independent standard normal variables. (The reader may wish to confirm this using Theorem 1 in Appendix A.) A drawback of the algorithm is that it involves calculation of expensive functions log, sin and cos. Also two random variables are generated even if only one is needed. From a standard normal variable, one can easily sample from the normal distribution N(µ, σ²); if z ∼ N (0, 1), then x = µ + zσ is the desired variable.

Random variables from other distributions such as the gamma and beta may be generated from special algorithms. These will not be discussed here; see Ripley (1987) and Ross (1997).

9.4 Generation of discrete random variables

(p. 302 ) 9.5 Simulating molecular evolution

9.5.1 Simulating sequences on a fixed tree

Here we consider generation of a nucleotide sequence alignment when the tree topology and branch lengths are given. The basic model assumes the same substitution process at all sites and along all branches, but we will also consider more complex models in which the evolutionary process may vary across sites or branches. We consider nucleotide models; amino acid or codon sequences can be generated using the same principles. Several approaches can be used and produce equivalent results.

9.5.1.1 Method 1. Sampling sites from the multinomial distribution of site patterns

If the substitution model assumes independent evolution at different sites in the sequence and all sites evolve according to the same model, data at different sites will have independent and identical distributions; they are said to be i.i.d. The sequence data set will then follow a multinomial distribution, with every site to be a sample point, and every site pattern to be a category (cell) of the multinomial. For a tree of s species, there are 4^s, 20^s, or 64^s possible site patterns for nucleotide, amino acid, or codon sequences, respectively. Calculation of the probability for every site pattern is explained in Section 4.2, which describes the calculation of the likelihood function under the model. A sequence alignment can thus be generated by sampling from this multinomial distribution. The result will be the numbers of sites having the site patterns, many of which may be zero if the sequence is short. If a pattern, say TTTC (for four species), is observed 50 times, one simply writes out 50 sites with the same data TTTC, either with or without randomizing the sites. Most phylogeny programs, especially those for likelihood and Bayesian calculations, collapse sites into patterns to save computation, since the probabilities of observing sites with the same data are identical. As this simulation method generates counts of site patterns, those counts should in theory be directly usable by phylogeny programs.

The multinomial sampling approach is not feasible for large trees because the number of categories becomes too large. However, it is very efficient for small trees with only four or five species, especially when combined with an efficient algorithm for sampling from the multinomial, such as the alias method.

9.5.1.2 Method 2. Evolving sequences along the tree

This approach ‘evolves’ sequences along the given tree, and is the algorithm used in programs such as Seq-Gen (Rambaut and Grassly 1997) and EVOLVER (Yang 1997a). First, we generate a sequence for the root of the tree, by sampling nucleotides according to their equilibrium distribution under the model: π_T, π_C, π_A, π_G. Every nucleotide is sampled independently from the discrete distribution (π_T, π_C, π_A, π_G). If the base frequencies are all equal, one can use the algorithm for discrete uniform distributions, which is more efficient. The sequence for the root is then allowed to (p. 303 ) evolve to produce sequences at the daughter nodes of the root. The procedure is repeated for every branch on the tree, generating the sequence at a node only after the sequence at its mother node has been generated. Sequences at the tips of the tree constitute the data, while sequences for ancestral nodes are discarded.

To simulate the evolution of a sequence along a branch of length t, calculate the transition probability matrix P(t) = {p_ij(t)} (see Sections 1.2 and 1.5), and then simulate nucleotide substitutions at every site independently. For example, if a site is occupied by C in the source sequence, the nucleotide in the target sequence will be a random draw from the discrete distribution of T, C, A, G, with probabilities p_CT(t), p_CC(t), p_CA(t), p_CG(t), respectively. This process is repeated to generate all sites in the target sequence. As the transition probabilities apply to all sites in the sequence for the branch, one has to calculate them only once for all sites for that branch.

9.5.1.3 Method 3. simulating the waiting times of a Markov chain

This is a variation to method 2. One generates a sequence for the root, and then simulates the evolution of any site along any branch as follows. Suppose the branch length is t, and the rate matrix of the Markov chain is Q = {q_ij}. Let q_i = −q_ii = Σ_j≠i q_ijbe the substitution rate of nucleotide i. Suppose the site is currently occupied by nucleotide i. Then the waiting time until the next substitution event has an exponential distribution with mean 1/q_i. We draw a random waiting time s from the exponential distribution. If s > t, no substitution occurs before the end of the branch so that the target sequence has nucleotide i at the site as well. Otherwise a substitution occurs and we decide which nucleotide the site changes into. Given that the site with nucleotide i has experienced a change, the probability that it changes into nucleotide j is q_ij/q_i, and we can sample j from this discrete distribution. The remaining time for the branch becomes t − s. We then draw a random waiting time from the exponential distribution with mean 1/q_j. The process is repeated until the time for the branch is exhausted.

This simulation procedure is based on the following characterization of a continuous-time Markov chain (Fig. 9.2). The waiting time until the next transition (change) is exponential with mean 1/q_i. If we ignore the waiting times between transitions, the sequence of states visited by the process constitutes a discrete-time Markov (p. 304 ) chain, which is called the jump chain. The transition matrix of the jump chain is given as (9.7) Note that every row sums to 1.

Fig. 9.2 Characterization of the Markov process as exponential waiting times and a jump chain. The waiting times until the next substitution s₁, s₂, and s₃ are independent random variables from the exponential distributions with means 1/q_T, 1/q_C, and 1/q_A, respectively, where q_T, q_C, and q_A are the substitution rates of nucleotides T, C, and A, respectively.

The algorithm of simulating exponential waiting times and the jump chain may be applied to the whole sequence instead of one site. The total rate of substitution is the sum of the rates across sites, and the waiting time until a substitution occurs at any site in the whole sequence has an exponential distribution with the mean equal to the reciprocal of the total rate. If a substitution occurs, it is assigned to sites with probabilities proportional to the rates at the sites.

An advantage of this simulation procedure is that it does not require calculation of the transition-probability matrix P(t) over branch length t, as both the waiting times and the transition matrix for the jump chain are fully specified by the instantaneous rates. As a result, this procedure can be adapted to simulate more complex sequence changes such as insertions and deletions. One simply calculates the total rate of all events (including substitutions, insertions, and deletions) for the whole sequence, and simulates the exponential waiting time until the next event. If an event occurs before the end of the branch, one assigns the event to a site and to an event type (a substitution, insertion, or deletion) with probabilities in proportion to the rates of those events at the sites.

9.5.1.4 Simulation under more complex models

The methods discussed above can be modified to simulate under more complex models, for example, to allow different substitution parameters among branches (such as different transition/transversion rate ratio κ, different base frequencies, or different ω ratios). The multinomial sampling approach (method 1) applies as long as the site pattern probabilities are calculated correctly under the model. The approaches of evolving sequences along branches on the tree, either by calculating the transition probabilities (method 2) or by simulating the waiting times and the jump chain (method 3) are also straightforward; one simply use the appropriate model and parameters for the branch when simulating the evolutionary process along that branch.

One may also simulate under models that allow heterogeneity among sites. We will consider as an example variable substitution rates but the approaches apply to other kinds of among-site heterogeneity. There are two kinds of models that incorporate rate variation among sites (see Subsections 4.3.1 and 4.3.2). The first is the so-called (p. 305 ) fixed-sites model, under which every site in the sequence belongs to a predetermined site partition. For example, one may simulate five genes evolving on the same tree but at different rates r₁, r₂, …, r₅. The transition probability matrix for gene k is p_ij(tr_k). Sites in different genes do not have the same distribution, but within each gene, the sites are i.i.d. Thus one can use either of the three methods discussed above to simulate the data for each gene separately and then merge them into one data set.

A second kind of heterogeneous-sites model are the so-called random-sites models. Examples include the gamma models of variable rates for sites (Yang 1993, 1994a), the codon models of variable ω ratios among sites (Nielsen and Yang 1998; Yang et al. 2000), and the covarion-like models (Galtier 2001; Huelsenbeck 2002; Guindon et al. 2004). The rates (or some other features of the substitution process) are assumed to be random variables drawn from a common statistical distribution, and we do not know a priori which sites have high rates and which have low rates. Data at different sites are i.i.d. The approach of multinomial sampling can be used directly, although the site pattern probabilities have to be calculated under the heterogeneous-sites model. One may also sample the rate for each site and apply the method of simulating evolution along branches. If a continuous distribution is used, one should in theory calculate the transition probability matrix P(t) for every site and every branch. If a few site classes are assumed (as in the discrete-gamma model), one may sample rates for sites first, and then simulate data for the different rate classes separately, perhaps followed by a randomization of the sites. Note that under the random-sites model, the number of sites in any site class varies among simulated replicates, whereas in the fixed-sites model, the number is fixed.

(p. 306 ) 9.6 Exercises

9.1 Write a small simulation program to study the birthday problem. Suppose that there are 365 days in a year and that one’s birthday falls on any day at random. Calculate the probability that at least two people out of a group of k = 30 people have the same birthday (that is, they were born on the same day and month but not necessarily in the same year). Use the following algorithm. (The answer is 0.706.)

1. 1. Generate k = 30 birthdays, by taking 30 random draws from1, 2, …, 365.

2. 2. Check whether any two birthdays are the same.

3. 3. Repeat the process 10⁶ times and calculate the proportion of times in which at least two out of 30 people have the same birthday.

9.2 Monte Carlo integration (Subsection 5.3.1). Write a small program to calculate the integral f (x) in the Bayesian estimation of sequence distance under the JC69 model, discussed in Subsection 5.1.2. The data are x = 90 differences out of n = 948 sites. Use the exponential prior with mean 0.2 for the sequence distance θ. Generate N = 10⁶ or 10⁸ random variables from the exponential prior: θ₁, θ₂, …, θ_n, and calculate (9.8) Note that the likelihood f (x|θ_i) may be too small to represent in the computer, so scaling may be needed. One way is as follows. Compute the maximum log likelihood ${\ell }_m=\log \left\{f(x|\widehat{\theta })\right\}$, where $\widehat{\theta }=0.1015$ is the MLE. Then multiply f (x|θ_i) in equation (9.8) by a big number e^−ℓm so that they are not all vanishingly small before summing them up; that is, (9.9)

9.3 Write a small simulation program to study the optimal sequence divergence when two sequences are compared to estimate the transition/transversion rate ratio κ under theK80 model. Assume κ = 2 and use a sequence length of 500 sites. Consider several sequence distances, say, d = 0.01, 0.02, …, 2. For each d, simulate 1000 replicate data sets under the K80 model and analyse it under the same model to estimate d and κ using equation (1.11). Calculate the mean and variance of the estimate $\widehat{k}$ across replicate data sets. Each data set consists of a pair of sequences, which can be generated using any of the three approaches discussed in Subsection 9.5.1.

9.4 Long-branch attraction by parsimony. Use the JC69 model to simulate data sets on a tree of four species (Fig. 9.3a), with two different branch lengths a = 0.1 and b = 0.5. Simulate 1000 replicate data sets. For each data set, count the sites with the three site patterns xxyy, xyxy, xyyx, and determine the most parsimonious tree. To simulate a data set, reroot the tree at an interior node, as in, say, Fig. 9.3(b). Generate (p. 307 ) a sequence for the root (node 0) by random sampling of the four nucleotides, and then evolve the sequence along the five branches of the tree. You may also use the approach of multinomial sampling. Consider a few sequence lengths, such as 100, 1000, and 10 000 sites.

Fig. 9.3 (a) A tree of four species, with three short branches (of length a) and two long branches (with length b) for simulating data to demonstrate long-branch attraction. (b) The same tree rerooted at an ancestral node for simulation.

9.5 A useful test of a new and complex likelihood program is to generate a few data sets of very long sequences under the model and then analyse them under the same model, to check whether the MLEs are close to the true values used in the simulation. As MLEs are consistent, they should approach the true values when the sample size (sequence length) becomes larger and larger. Use the program written for Exercise 9.4 to generate a few data sets of 10⁶, 10⁷, or 10⁸ sites and analyse them using a likelihood program (such as PHYLIP, PAUP, or PAML) under the same JC69 model, to see whether the MLEs of branch lengths are close to the true values. Beware that some programs may demand a lot of resources to process large data sets; save your important work before this exercise in case of a computer crash.