Monte Carlo Evaluation of the Likelihood for Ne From Temporally Spaced Samples

FORMULATION OF THE MODEL AND MONTE CARLO

The model: The data are genetic samples collected at different generations. The first sample is collected at generation 0 and the last sample at generation T. Any samples drawn at intervening generations may be evenly or irregularly spaced in time. For notational simplicity, we assume for now that individuals are genotyped at a single locus, though we describe later the extension to multiple, independently segregating loci. The data include K different allelic types, indexed by k = 1,..., K. The allele frequencies observed in samples taken from different generations will differ due to genetic drift and sampling variation.

Let Y_t = (Y_t,1,..., Y_t,K) be the counts of the K different allelic types in the sample at generation t, and let S_t denote the number of diploid individuals in the sample. We assume that the samples were taken from a Wright-Fisher population of size N_e and denote the unobserved population allele counts at generation t by X_t = (X_t,1,..., X_t,K), with Σk=1KXt,k=2Ne . By the formulation of the Wright-Fisher model, the X_t form a Markov chain in time, with transitions defined by multinomial probabilities depending on N_e, PNe(Xt∣X0,…,Xt−1)=PNe(Xt∣Xt−1)=(2Ne)!∏k=1K[Xt−1,k∕(2Ne)]Xt,kXt,k!,t=1,2,…. (1)

The genetic sample at a time t is assumed to be drawn with replacement from the copies of alleles present in the population at time t. This is equivalent to drawing the sample Y_t from a very large gamete pool produced by the population at time t: sampling plan II of Waples (1989). This type of sampling applies to many organisms, especially those species with high fecundity that may be sampled as juveniles, or those that may be sampled (preferably noninvasively) as adults in populations having census sizes considerably larger than their effective sizes (Waples 1989). The sample allele counts Y_t, given the latent variable X_t, are conditionally independent of all the other variables and follow the multinomial distribution depending on the parameter N_e, the sample size S_t, and X_t, PNe(Yt∣Xt)=(2St)!∏k=1K[Xt,k∕(2Ne)]Yt,kYt,k!, (2) when S_t > 0. If there is no sample taken from the population at generation t, then S_t ≡ 0, and we define P_Ne(Y_t|X_t) ≡ 1.

Such a system forms a hidden Markov chain with the dependence structure shown in the directed graph of Figure 1. The allele counts in the population when the first sample is drawn, X₀, are nuisance parameters. To avoid estimating X₀ and the associated consistency problems with the maximum-likelihood estimator (Neyman and Scott 1948), we consider the integrated likelihood, assuming a uniform prior distribution, π(X₀), on the population allele counts at time 0. The likelihood for N_e is the probability of the data Y = (Y₀,..., Y_T) given the parameter N_e. The probability of Y is the sum of the joint probability of Y and the latent variables X = (X₀,..., X_T) over the space of all X: PNe(Y)=∑XPNe(Y,X)=∑X0,…,XT(π(X0)PNe(Y0∣X0)∏t=1TPNe(Xt∣Xt−1)PNe(Yt∣Xt)). (3) For the case of K = 2 and N_e small, the likelihood in (3) may be computed exactly. Williamson and Slatkin (1999) effected the summation in (3) in terms of multiplication of transition probability matrices. The dimension of the square matrices is (N_e - 1)!/[(N_e - K)!(K - 1)!], which increases rapidly with N_e and K. We note that the hidden Markov form of the system allows a more efficient computation of the likelihood using the algorithm of Baum (1972). Nonetheless, exact evaluation for multiple alleles would still require prohibitively large amounts of computation and storage. An alternative is to estimate P_Ne(Y) by Monte Carlo.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

—A directed graph showing the dependence structure of the components of X and Y. The Y_T’s are observations of a hidden Markov chain. The graph shown represents a situation where samples were taken at generations 0, 2, t, and T, and no samples were taken at generations 1 and t - 1.

Monte Carlo evaluation: For likelihood inference, we must evaluate P_Ne(Y) for a number of different values of N_e. Expressing this probability as an expectation with respect to the distribution of X gives PNe(Y)=∑XPNe(Y,X)=∑XPNe(Y∣X)PNe(X)=ENe(PNe(Y∣X)). (4) In this form the expectation would be taken over the marginal probabilities of X, and it could be estimated by Monte Carlo as PNe(Y)≈1m∑i=1mPNe(Y∣X(i)) (5) for large m, with X⁽ⁱ⁾ being the ith realization from the marginal distribution of X. Such a naive scheme fails, however, because P_Ne(Y|X⁽ⁱ⁾) varies greatly over the values of X realized from their marginal distribution, resulting in enormous Monte Carlo variance.

Instead, we pursue a more efficient Monte Carlo approximation by using importance sampling (Hammersley and Handscomb 1964). We express P_Ne(Y) as an expectation with respect to a different distribution PNe∗(X) having the property that PNe∗(X)>0 0 for all X such that P_Ne(Y, X) > 0. Thus, we have PNe(Y)=∑XPNe(Y,X)PNe∗(X)PNe∗(X)=ENe∗(PNe(Y,X)PNe∗(X)), (6) where ENe∗ indicates that the expectation is over the space of X weighted by the distribution PNe∗(X) . The expectation (6) may be estimated by Monte Carlo, giving PNe(Y)≈P∼Ne(Y)=1m∑i=1mPNe(Y,X(i))PNe∗(X(i)) (7) for large m, where X⁽ⁱ⁾ is the ith realization of X drawn from PNe∗(X) . The Monte Carlo variance of P˜_Ne(Y) is made small when PNe(Y,X)∕PNe∗(X) varies little across the possible values of X and would be minimized if PNe∗(X) were exactly proportional to P_Ne(Y, X). Such a distribution of X would, by definition, be the conditional distribution P_Ne(X|Y). Unfortunately, for the same reasons that P_Ne(Y) cannot be computed exactly, it is infeasible to compute P_Ne(X|Y). Nonetheless, the Monte Carlo variance of P˜_Ne(Y) will be reduced to the extent that PNe∗(X) resembles P_Ne(X|Y). We now describe a method for rapid simulation of X⁽ⁱ⁾’s from a distribution PNe∗(X) that is close to P_Ne(X|Y). As is required for the importance sampling, it is also possible to compute PNe∗(X(i)) quickly for each X ⁱ⁾ generated.

Sampling from PNe∗(X) by a forward-backward method: Baum et al. (1970) describe a method applicable to general, hidden Markov chains for realizing latent variables, such as X = (X₀,..., X_T), from their exact conditional distribution given the observed variables, such as Y = (Y₀,..., Y_T). Their algorithm first employs a “forward step” in which the conditional probability distributions of each X_T, given the observed variables up to and including Y_T, are recursively computed and stored using the relation P(Xt∣Y0,…,Yt)∝∑Xt−1P(Xt−1∣Y0,…,Yt−1)P(Xt∣Xt−1)P(Yt∣Xt), (8) which is normalized by the sum of that quantity over all the values of X_t. The last such conditional distribution computed is P(X_T|Y₀,..., Y_T). The “backward step” begins with simulating a value XT(i) from this distribution (where, as before, the superscript (i) indicates a realized value of a random variable). One then proceeds backward, realizing XT−1(i) from its conditional distribution given all of the observed variables, Y, and XT(i) . In similar fashion, one realizes XT−2(i) and so forth back to X0(i) . In this backward phase, each Xt(i) is simulated from its conditional distribution given all the data Y and all of the components of X that have been realized so far. That is, Xt(i) is drawn from P(Xt∣Y0,…,YT,Xt+1(i),…,XT(i)). (9) Because of the conditional independence structure in a hidden Markov chain, (9) reduces to P(Xt∣Y0,…,Yt,Xt+1(i)) , which may be computed using the distributions stored during the forward step by the relation P(Xt∣Y0,…,Yt,Xt+1(i))∝P(Xt∣Y0,…,Yt)P(Xt+1(i)∣Xt). (10) At the end of the backward step, it is thus clear that the resulting realization (X0(i),…,XT(i)) is from the conditional distribution of X given Y.

An approximation for multiple alleles: In our application, with multiple alleles at a locus, since there are so many possible states that each X_t may take, the above procedure is computationally infeasible. However, we make use of the Baum et al. (1970) algorithm in spirit, employing two alterations to make it feasible to simulate from PNe∗(X) and to compute its value. We emphasize that although the method described below involves a series of “approximations” by which PNe∗(X) , differs from P_Ne(X|Y), the final sampling and computation of PNe∗(X) is exactly from the PNe∗(X) as constructed, so its use in (7) gives a true Monte Carlo estimate.

The first approximation is to perform the forward-backward cycle separately for each allele. To describe this, we introduce some more notation. Denote by X_(k) the vector (X_0,k,..., X_T,k) of latent counts of the kth allele from time t = 0 to t = T. Similarly we define Y_(k) = (Y_0,k,..., Y_T,k). To do the forward-backward cycle separately over alleles we first focus on allele 1, simulating X(1)(i) by the forward-backward mechanism as if the data were on a diallelic locus with observed counts Y₍₁₎ from samples of size S₀,..., S_T through time. Once we have realized X(1)(i) we update the sizes of the population and the sample. Thus we define the updated population size vector 2N(2)∗=(2Ne−X0,1(i),…,2Ne−XT,1(i)) and an updated sample size vector 2S(2)∗=(2S0,2∗,…,2ST,2∗)=(2S0−Y0,1,…,2ST−YT,1) , in effect removing the first allelic type from the remainder of the data and the population. We then use the forward-backward mechanism again to simulate X(2)(i) as though the data were counts Y₍₂₎ from a diallelic locus drawn from a population with sizes that change over time N(2)∗ and sample sizes S(2)∗ . This continues sequentially over alleles, updating population sizes and sample sizes as above: 2N(k)∗←(2N(k−1)∗−X(k−1)(i)) and 2S(k)∗←(2S(k−1)∗−Y(k−1)(i)) , until X_(K-1) has been realized, which also determines that X(K)←(2N(K−1)∗−X(K−1)(i)) . (Here and later we use the notation A ← B to mean “the value B is assigned to the variable A.”) At the end one has obtained a realized value X⁽ⁱ⁾, which may be used in (7).

PNe∗(X) using a continuous approximation: Although realizing alleles sequentially, as above, greatly reduces the number of terms required to use (8) and (10), the method would still involve a prohibitive amount of summation over binomial probabilities. Thus we construct PNe∗(X) , employing a normal approximation to binomial probabilities, which replaces all such sums by analytically tractable integrals. Recall that if W ∼ Binomial(n, p), then the transformed variable sin^-1(W/n)^1/2 is approximately normally distributed with variance 1/(4n). Note that this quantity does not depend on p. Hence we use this transformation to define the quantities ϕt,k=sin−1[Yt,k∕(2St,k∗)]1∕2 when St,k∗>0 , and θt,k=sin−1[Xt,k∕2Nt,k∗)]1∕2 . By realizing the continuous values θt,k(i) in a forward-backward framework within a continuous setting, the computational demands are greatly reduced. And then, by transforming each θt,k(i) back into the appropriate discrete Xt,k(i) we have a way to realize X⁽ⁱ⁾ from PNe∗(X) and to compute the probability PNe∗(X(i)) . The details of this procedure are given in the appendix. We use it to compute the Monte Carlo estimate P˜_Ne(Y) using (7).

Monte Carlo variance and multiple loci: The quantity P∼Ne(Y) is only an estimate of the true value P_Ne(Y). By the central limit theorem, for large m, P_Ne(Y) will be approximately normally distributed (Hammersley and Handscomb 1964) with mean P˜_Ne(Y) and a variance that may be approximated without bias by the quantity Var^(P∼Ne(Y))=1m(m−1)∑i=1m(PNe(Y,X(i))PNe∗(X(i))−P∼Ne(Y))2. (11) These facts may be used to obtain a confidence interval estimate around each P˜_Ne(Y).

The ability to estimate N_e typically requires data from many loci. The extension to data on J independently segregating loci, indexed by j = 1,..., J, is straightforward—each locus is treated separately, and the estimated likelihoods from each locus are multiplied together. Thus, let P˜_Ne,_j(Y) be the Monte Carlo likelihood estimate from the data on the jth locus. The Monte Carlo likelihood estimate using all the loci is then P∼NeJ(Y)=∏j=1JP∼Ne,j(Y). (12) This requires that the initial allele counts have independent prior distributions, π(X₀). An implicit assumption of (12) is consequently that the loci used are in linkage equilibrium at t = 0. P∼NeJ(Y) will also have an approximately normal distribution. An unbiased estimator for its Monte Carlo variance (derived in Equation A9 in the appendix) is Var^(P∼NeJ(Y))=∏j=1J(P∼Ne,j(Y))2−∏j=1J([P∼Ne,j(Y)]2−Var^(P∼Ne,j(Y))). (13) This can be used to compute a confidence interval estimate around P∼NeJ(Y) .

When displaying the Monte Carlo likelihood curve it is preferable to plot the log-likelihood values, logP∼NeJ(Y) , for different values of N_e. In this case, the endpoints of the confidence intervals may be similarly log-transformed.

SIMULATED AND REAL DATASETS

We demonstrate the method by computing log-likelihood curves for N_e from three different datasets. First, to verify that the method gives correct results we apply it to a simple simulated dataset (dataset 1) for which it is possible to compute the likelihood exactly. We simulated samples of 20 diallelic loci from 100 diploid individuals at generations 0, 6, and 12 from a Wright-Fisher population of 25 diploid individuals. For each locus, the initial allele frequency in the population at time zero was an independently drawn uniform real number between 0 and 1. The log-likelihood for N_e given these data was estimated for values between 10 and 52, in steps of 2, using m = 20,000 realizations of X from PNe∗(X) for each locus and each N_e.

We simulated a second dataset (dataset 2) to see how the method performed with multiallelic markers taken from a Wright-Fisher population. The dataset included three samples of 100 diploids for 12 five-allele loci at generations 0, 4, and 8 from a population of 50 diploids. The allele frequencies at each locus in generation 0 for these simulations were independently drawn from a uniform Dirichlet density with five components. For this dataset, the log-likelihood was computed for values of N_e between 20 and 100 in increments of 4 using m = 50,000 realizations of X for each locus and each value of N_e.

Finally, we computed a log-likelihood curve for N_e given data on a population of Drosophila in Begon et al. (1980). These data were analyzed using F-statistics by Begon et al. (1980) as well as by Pollak (1983). They observed allele frequencies in three samples at each of nine enzyme loci. The first two samples were taken a little more than 1 yr apart, and the third sample was taken some 8 mo later. Though the natural populations do not have discrete generations, they have been modeled previously by Begon et al. and Pollak as populations with discrete generations. Because of the different growth rates of flies during different seasons, seven generations separate the first two samples, while only two generations separate the second two samples (Begonet al. 1980). The sample sizes for all loci were the same, with larger sample sizes taken in the latter sampling periods. The sample sizes were S₀ = 190, S₇ = 250, and S₉ = 335 flies. Pollak (1983) notes that since Begon et al. (1980) sampled adult flies, their sampling scheme is closer to what is known in the literature as sampling scheme I than it is to sampling scheme II. However, as discussed by Waples (1989), the probability models underlying the two different sampling schemes are very similar when the actual size of the population is much larger than the effective size of the population. This is the case with these Drosophila. Begon et al. (1980) report census sizes in the tens of thousands of flies, while the estimated N_e is orders of magnitude smaller. Because of this, it is still reasonable to analyze the data using the likelihood method we have developed here.

The data appear as allele frequencies in Table 1 of Begon et al. (1980). Unfortunately the allele frequencies at the Pgm locus are misreported there and fail to sum to one. We thus used only the remaining eight loci. Of these eight, three had three alleles, two had four alleles, two had five alleles, and one had six alleles. We evaluated P∼NeJ(Y) at values of N_e between 200 and 1200 in increments of 50, with two more points (N_e = 425 and N_e = 475) included near the peak of the likelihood curve. For each point we used m = 500,000 realizations of X.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

—Log-likelihood curves estimated by Monte Carlo from datasets 1 and 2. The values of N_e at which the likelihood was computed are indicated by vertical lines above the horizontal axis in each figure. The log-likelihood values are connected by a solid line. Vertical bars intersecting the solid line indicate 90% confidence intervals around log computed using the Monte Carlo variance estimate (13). The endpoints of the confidence intervals are connected by dashed lines. (a) Dataset 1 is simulated data from 20 diallelic loci. (b) Dataset 2 is simulated data from 12 loci with five alleles each.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

—Log-likelihood curves from the data of Begon et al. (1980) estimated by Monte Carlo. The format of the plot is as for Figure 2.

APPENDIX

Using θ_t,k and ϕ_t,k in a continuous setting: We define the random variables ϕt,k=sin−1[Yt,k∕(2St,k∗)]1∕2 when St,k∗>0 , and θt,k=sin−1[Xt,k∕(2Nt,k∗)1∕2 when Nt,k∗>0 . These quantities have an approximate normal distribution, which is independent of their means. We use them in our construction of the importance-sampling function P_Ne(X). Below, we concentrate on their use for realizing X(k)(i) , keeping in mind that if k > 1 then we will have already realized X(k−1)(i) , and we use the updated population and sample sizes N(k)∗ and S(k)∗ . If k = 1 then N(1)∗=(Ne,…,Ne) and S(1)∗=(S0,…,ST) , respectively.

The forward step: Following Cavalli-Sforza and Edwards (1967), if θ_t-1,k is normally (+) distributed with mean μ_t-1 and variance σt−12 , then, after a generation of genetic drift in a population of Nt,k∗ diploids, θ_t,k has an approximate normal distributions with mean μ_t-1 and variance σt2=σt−12+1∕(8Nt,k∗) . If there are data Y_t,k from a sample of size S_t,k at time t, then ϕ_t,k has an approximate normal distribution with mean θ_t,k and variance 1∕(8St,k∗) , so, given that θt,k∼N(μt,σt2) , the conditional distribution of θ_t,k given ϕ_t,k is also normal. These relations form the basis of a continuous approximation for doing the forward step. For the purpose of realizing X we assume that the uniform prior on X₀ is equivalent to a diffuse prior on θ_0,k. Therefore θ0,k∣ϕ0,k∼N(μ0,σ02) with μ₀ = ϕ_0,k σ02=1∕(8S0,k∗) . With that as a starting point, we work iteratively forward in time, assigning values μt←μt−1 (A1) σt2←σt−12+1∕(8Nt,k∗) (A2) if St,k∗=0 . If St,k∗>0 , however, then one first computes μ_t and σ² as in (A1) and (A2), but then further updates the values to reflect the information in the sample at time t: μt←μt∕(8St,k∗)+σt2ϕt,k1∕(8St,k∗)+σt2 (A3) σt2←σt2∕(8St,k∗)1∕(8St,k∗)+σt2. (A4) This is analogous to computing a posterior distribution from a normal prior and normal data (see, for example, Gelmanet al. 1996, p. 43).

Carrying this out until t = T gives values for the mean and variance of θ_T,k given ϕ_0,k,..., ϕ_T,k, assuming they follow a normal distribution. In fact, for each t, it gives us the parameters for the normal distribution of θ_t,k conditional on ϕ_r,k for all r ≤ t. We are thus in a position to realize θt,k(i)’s in the backward step and transform those θt,k(i) ’s back into the Xt,k(i) ’s that we need.

The backward step: The backward step is more complicated than the forward step, because after realizing each value of θt,k(i) we must transform it into the discrete value Xt,k(i) that we require. This transformation process requires some extra bookkeeping to ensure that we do not waste time realizing X⁽ⁱ⁾’s that are incompatible with the data. This is described in the next section of the appendix. We first realize the value θT,k(i) from a N(μT,σT2) distribution. Then we transform that to the realization XT,k(i) by a many-to-one map M , which has two effects: the first is that of folding and translating the distribution of θT,k so that it is bounded between 0 and π/2, mapping θT,k(i)∈(−∞,∞) to a value θT,k∗∈[0,π∕2] . The second is transforming that θT,k∗ into the appropriate value XT,k(i) (see the next section).

Working backward, each θT,k(i) for t = T - 1 down to t = 0, is realized from a N(μt,σt2) distribution and then transformed into the corresponding θT,k∗ and Xt,k(i) by M . In keeping with (10), before θt,k(i) is realized, μ_t and σt2 must be appropriately updated, on the basis of the values of μ_t and σt2 stored during the forward step and the realized value θt+1,k(i) . This involves making the assignments μt←μt∕(8Nt+1∗)+σt2θt+1,k∗1∕(8Nt+1∗)+σt2 (A5) σt2←σt2∕(8Nt+1∗)1∕(8Nt+1∗)+σt2 (A6) in the order as written.

Computing the probability PNe∗(X(i)) : By carrying out the forward-backward steps above on the first allele, the realization X(1)(i) is obtained. Then, N(2)∗ and S(2)∗ are computed and used in the forward-backward steps to obtain X(2)(i) . Executing these steps for all the alleles yields the realization X⁽ⁱ⁾, which is used in (7). P_Ne(Y, X⁽ⁱ⁾) in (7) is easily computed using the expansion shown between the large parentheses in (3).

It remains only to compute PNe∗(X(i)) , which can be done by recording the probability of realizing each component Xt,k(i) . Although this probability depends on the values of μ_t, σt2 , N(k)∗ , and several bookkeeping variables, we denote it here simply by P(Xt,k(i)) . (The actual function P is described later in this appendix.) As long as the realization of X(k)(i) over alleles occurs in the same order over k (k = 1, 2,..., K) for each i, then PNe∗(X(i))=∏k=1K∏t=1TP(Xt,k(i)). (A7)

Details of M : The fact that we are realizing X(k)(i) ’s one allele at a time requires that we do some extra bookkeeping to keep our importance-sampling scheme efficient. Primarily, we must avoid realizing X⁽ⁱ⁾’s for which P_Ne (Y, X⁽ⁱ⁾) = 0. Potential problems arise because by the method we use to realize values from PNe∗(X),Xt,k may only take values between 0 and 2Nt,k∗ , inclusive. If 2Nt,k∗=0 at any value of t, then for any s>t,Xs,k(i) must also be 0. To avoid situations in which this leads to P_Ne (Y, X⁽ⁱ⁾) being 0 (as when Xt,k(i)=0 and Y_t,k > 0) we introduce the following scheme and additional notation: δt,k={1ifXi,k(i)=0impliesPNe(Y,X(i))=00otherwise}γt,k=minr<t2Nr,k∗κt,k=the number of allelic subscriptsℓ:k<ℓ≤Ksuch thatYr,ℓ>0for at least oner≥t. (A8) Knowing the above quantities, we can define the function M . In the remainder of this section and in the following one we drop the t and k subscripts for clarity.

With N^* and γ positive integers, δ ∊ {0, 1}, and κ ∊ {0, 1,..., min(2N^* - δ, γ - δ}, let M(θ;N∗,δ,γ,κ):R1→{δ,…,2N∗−κ}×[0,π∕2] be the many-to-one map that takes a realization of θ ∊ (-∞, ∞) to the ordered pair (X, θ^*), where X is an integer such that δ ≤ X ≤ 2N^* - κ, and θ^* is a real number between 0 and π/2, inclusive. M may be described by the following pseudocode. We first define the quantities L = sin^-1(0.5/(2N^*))^1/2 and H={sin−1[(2N∗−κ+0.5)∕(2N∗)]1∕2,κ≥1sin−1[(2N∗−0.5)∕(2N∗)]1∕2,κ=0.} Then,

if (δ = 2N^* - κ or δ = γ - κ = 0) then θ^* ← 0

else if (L < θ < H) then θ^* ← θ

else if (θ < L)

and if (δ = 0) then θ^* ← θ

else if (δ = 1) then θ_[L] ← 2L - θ (this is reflection around θ = L), and then

if (L ≤ θ_[L] < H) then θ^* ← θ_[L]

else we know θ_[L] ≥ H, and we consider the sequence θ_[i] = i(L - H) + θ_[L], i = 1, 2,..., and we assign θ^* ← θ_[i^*], where i^* is the least i such that L < θ_[i] < H. (The sequence θ_[i] represents successive translation leftward.)

else if (θ > H)

and if (κ = 0) then θ^* ← π/2

else if (κ > 1) then θ_[H] ← 2H - θ (this is reflection around θ = H), and then

if (L < θ_[H] < H) then θ^* ← θ_[H]

else we know θ_[H] < L and we consider the sequence θ_[j] = j(H - L) + θ_[H], j = 1, 2,..., and we assign θ^* ← θ_[j^*], where j^* is the least j such that L ≤ θ_[j] < H. (The sequence θ_[j] represents successive translation rightward.)

finally we use θ^*, making the assignment X←⌊2N∗sin2θ∗+0.5⌋ ,

where ⌊x⌋ denotes the largest integer ≤x. The reflections and translations are depicted graphically in Figure A1a.

• Download figure
• Open in new tab
• Download powerpoint

Figure A1.

—Figures representing and for 2N^* = 20. The normal curve is the density for θ. (a) Reflections and translations as described in the appendix. Long-dashed lines represent the curve after reflection through L or H, while the short-dashed lines represent the reflected curve after one or more successive translations. (b) If δ = 1 and κ = 1, then X⁽ⁱ⁾ is constrained to be in {1..., 2N^* - 1}. The shaded regions correspond to those values of θ for which X⁽ⁱ⁾ = 13 by . The total shaded area is equal to . (c) If δ = 0 then X⁽ⁱ⁾ may take the value 0. The shaded area shows . (d) If κ = 0 then X may take the value 2N^*. The shaded area shows ; 10, δ, γ, 0).

The probability Pμ,σ2(X=X(i);N∗,δ,γ,κ) of realizing X=X(i) : If θ∼N(μ,σ2) , and (X,θ∗)=M(θ;N∗,δ,γ,κ) , then we denote by Pμ,σ2(X=X(i);N∗,δ,γ,κ) the marginal probability that X = X⁽ⁱ⁾. The value of Pμ,σ2(X=X(i);N∗,δ,γ,κ) can be expressed using the notation from the above section. First, Pμ,σ2(X=X(i);N∗,δ,γ,κ)=0 if X⁽ⁱ⁾ < δ or X⁽ⁱ⁾ > 2N^* - κ, though such values of X⁽ⁱ⁾ should never occur from M anyway. Second, there are cases when M constrains X⁽ⁱ⁾ to be either 0 or 1 with probability 1. Hence if 2N^* - κ = δ or γ - κ = δ = 0 then Pμ,σ2(X=δ;N∗,δ,γ,κ)=1 .

If, on the other hand, 2N^* - κ > δ and γ - κ > 0, then for X⁽ⁱ⁾ = 0 and X⁽ⁱ⁾ = 2N^* we have Pμ,σ2(X=0;N∗,0,γ,κ)=P(−∞<θ<L)Pμ,σ2(X=2N∗;N∗,δ,γ,0)=P(H≤θ<∞), while for 0 < X⁽ⁱ⁾ < 2N^* - κ we define a = sin^-1[(X⁽ⁱ⁾ - 0.5)/2N^*)]^1/2 and b = sin^-1[(X⁽ⁱ⁾ + 0.5)/(2N^*)]^1/2 and have Pμ,σ2(X=X(i);N∗,δ,γ,κ)=P(a≤θ<b)+I{δ=1}P(a≤θ[L]<b)+I{κ>0}P(a≤θ[H]<b)+I{δ=1}∑i=1∞P(a≤θ[i]<b)+I{κ>0}∑j=1∞P(a≤θ[j]<b), where I{·} is the indicator function and P(a ≥ θ < b) is the probability that a N(μ,σ2) random variable is between a and b, namely ∫ab(2πσ2)−1∕2exp{[−(θ−μ)2]∕(2σ2)}dθ . We compute this probability by numerical integration in our programs. In practice, the infinite sums are approximated by summing the first several terms of the series, until the contribution of the next term is very small (e.g., <10^-7). Values of P for different values of δ and κ appear as shaded regions in Figure A1, b-d.

This folding and translating might seem to be a very involved process, but it is computationally much faster than realizing θ from a truncated normal distribution and computing the probability of X⁽ⁱ⁾ when θ is from such a distribution.

Multilocus Monte Carlo variance calculation: We derive an expression for the variance of a product of J independent random variables, W_j, j = 1,... J: Var(∏Wj)=E([∏Wj]2)−[E(∏Wj)]2(definition of variance)=E([∏Wj2])−[E(∏Wj)]2(powers distribute over products)=∏E(Wj2)−∏[E(Wj)]2(independence of theWj)=∏E(Wj2)−∏[E(Wj2)−Var(Wj)](definition of variance). Denoting P˜_Ne,j(Y) in (13), by W_j and taking the expectation gives the same result, verifying that the expression in (13) is unbiased for Var(P∼NeJ(Y) ).