Introduction
============

Methods for inferring the selection coefficient at a single genetic locus from time series data have been employed extensively in evolutionary studies of simple traits. Such methods track the frequency of an allele or Mendelian trait over multiple generations and infer the selection coefficient that best explains the observed frequency changes. Studies of selective pressures conducted using time series approaches have provided evidence for selective forces in natural populations and have helped to characterize the ways in which environmental factors influence evolution through selection ([@msw173-B18]; [@msw173-B6]; [@msw173-B53]; [@msw173-B35]; [@msw173-B50]; [@msw173-B22]; [@msw173-B24]; [@msw173-B8]; [@msw173-B26]; [@msw173-B3]; [@msw173-B42]; [@msw173-B9]; [@msw173-B32]).

Because random fluctuations in allele frequencies due to genetic drift are often small compared with changes due to selective pressures, it is common practice for studies to assume that allele frequencies change deterministically over time as they would in a population of infinite size according to well-known deterministic formulas of ([@msw173-B17], p. 424) and ([@msw173-B23], p. 840) or related expressions ([@msw173-B21]; [@msw173-B25]). However, because allele frequency trajectories can be heavily influenced by genetic drift when population sizes or selection coefficients are small, many methods have been developed to account for drift by explicitly modeling finite population sizes when inferring selection coefficients from observed allele frequency trajectories ([@msw173-B37]; [@msw173-B40]; [@msw173-B2]; [@msw173-B36]; [@msw173-B38]; [@msw173-B33]; [@msw173-B48]; [@msw173-B19]; [@msw173-B16]; [@msw173-B45]) and when testing hypotheses about selection versus drift ([@msw173-B18]; [@msw173-B57]; [@msw173-B55]; [@msw173-B39]; [@msw173-B14]).

Although it is commonly assumed that estimates of selection coefficients are likely to be improved by accounting for population size histories, the expected amount of improvement is not well characterized. Even in relatively small populations, allele frequencies and other evolutionary processes behave almost deterministically if the selection coefficient or allele frequency is sufficiently high ([@msw173-B43]), suggesting that methods that ignore drift might perform well under these conditions. Conversely, if drift is strong allele frequency trajectories can be noisy and the accuracy of methods that ignore drift may be comparable to that of methods that account for population size, as all methods are likely to perform poorly under these conditions ([@msw173-B20]).

If computationally fast methods that ignore drift are accurate, they could dramatically reduce the time required to infer selection coefficients in data sets with many loci. In addition to their computational efficiency, methods that ignore drift by assuming that the population size is infinite do not require estimates of effective population sizes, which can be difficult to obtain accurately. Moreover, ignoring drift can lead to simple formulas and inference procedures under complicated scenarios involving multiple populations with migration or multiple loci with recombination ([@msw173-B28]). Therefore, in light of the beneficial properties of methods that ignore drift and assume deterministic allele frequency trajectories, it is of interest to compare their accuracy to that of methods that account for population size histories.

The theoretical accuracy of methods for inferring selection coefficients can be difficult to derive analytically. Thus, to explore differences between methods that ignore or account for drift, one can take the approach of empirically comparing inferences made by estimators that either account for the true population size history or ignore the size history by assuming that populations are large and drift is negligible. This is the approach we take here. For our analyses, we consider maximum likelihood estimators of selection coefficients because they are typically quite accurate and have desirable statistical properties. Moreover, the majority of recently developed methods for inferring selection coefficients from time series data are maximum likelihood estimators, making them an important category of methods to evaluate.

To draw conclusions about the accuracy of maximum likelihood estimators, it is important to consider estimators based on exact likelihoods rather than approximations, so that differences in estimates can be attributed entirely to whether a method ignores or accounts for drift. Although several approximate approaches have been developed for computing the likelihood of a selection model given time series allele frequency data, only three existing methods compute probabilities that are exact under a widely accepted model. In particular, the methods of [@msw173-B2] and [@msw173-B48] compute exact probabilities under the diffusion approximation of the Wright--Fisher process. However, these methods do not model time-varying population size histories. The third inference method based on an exact likelihood considers time-varying population size histories under the diffusion approximation of the Wright--Fisher process ([@msw173-B45]); however, it uses an MCMC algorithm to perform Bayesian inference that is not easily incorporated into a unified inference algorithm that allows us to directly compare inferences made by estimators that model the true population history with those that assume a population of infinite size without drift. No existing method computes the exact probability of an allele frequency trajectory under the discrete Wright--Fisher model, as the matrix powers required for such a method are considered to be computationally inefficient.

Here, we derive the exact probability of an allele frequency trajectory in a population of piecewise constant size under two classical models: the discrete Wright--Fisher model and the diffusion approximation of the Wright--Fisher process. We then use maximum likelihood estimators obtained using these probabilities to explore how ignoring or accounting for the true population history affects estimates of selection coefficients.

Our results are useful for understanding when point estimates obtained using estimators based on deterministic allele frequency trajectories are likely to be accurate and when accounting for the true population history could improve these estimates. Our results have implications for the interpretation of existing estimates of selection coefficients and for the use of deterministic estimators in future studies. The results can also help guide the development of demography-aware estimators of selection coefficients by identifying scenarios under which such estimators are likely to improve inference accuracy.

Results
=======

To compare the performance of estimators that ignore or account for drift, we inferred selection coefficients from allele frequency trajectories simulated under a variety of population histories of time-varying size.

The Population Model
--------------------

In all of our analyses, we considered a single biallelic locus with alleles labeled *a* and *A* evolving under selection and recurrent mutation in a panmictic population comprised of *L* different epochs $\ell = 1,\ldots,L$, each with constant size $N_{\ell}$ diploid individuals ([fig. 1](#msw173-F1){ref-type="fig"}). Epoch ℓ corresponds to the time interval $\lbrack\tau_{\ell - 1},\tau_{\ell}\rbrack$, where time is measured continuously in units of generations and we define $\tau_{0} \equiv 0$. By varying the population sizes $N_{\ell}$ across epochs, it is possible to model a variety of size-change patterns including exponential growth, bottlenecks, and rapidly oscillating population sizes. Fig. 1.Diagram of the model. An allele at a single locus evolves in a population of piecewise constant size with $L = 5$ epochs spanning the time periods $\lbrack\tau_{0},\tau_{1}\rbrack,\ldots,\lbrack\tau_{L - 1},\tau_{L}\rbrack$, where $\tau_{0} \equiv 0$. Samples of sizes $n_{1},\ldots,n_{k}$ haplotypes are taken at times $t_{1},\ldots,t_{K}$.

Within epoch $\ell$, all mutation and selection parameters are assumed to be constant. In particular, we assume that the per-generation probability that allele *a* mutates to allele *A* is $u_{aA}^{(\ell)}$ and the per-generation probability that allele *A* to *a* is $u_{Aa}^{(\ell)}$. The three possible genotypes, *aa*, *aA*, and *AA*, have relative fitnesses given by $w_{AA}^{(\ell)} = 1 + s_{\ell},\, w_{aA}^{(\ell)} = 1 + h_{\ell}s_{\ell}$, and $w_{aa}^{(\ell)} = 1$ in epoch $\ell$, where $s_{\ell}$ is the selection coefficient and $h_{\ell}$ is the dominance parameter.

We denote the collection of model parameters in epoch $\ell$ by $\Theta_{\ell}$ and the set of parameters across all epochs by Θ. It will also be convenient to denote the value of the model parameters at time *t* by *N~t~*, $u_{aA}^{(t)},\, u_{Aa}^{(t)}$, *s~t~*, and *h~t~*, where *t* is measured continuously in units of generations. The epoch in which time *t* lies will be denoted by $\ell_{t}$ and the epoch in which sampling event *k* lies will be denoted by $\ell_{k}$. It will be clear from the context whether the subscript on $\ell$ refers to a time or a sampling event.

We denote the population-wide number of copies of allele *A* in generation *t* by *c~t~* and the population-wide frequency of allele *A* by *y~t~*. In practice, we do not observe the true population count of allele *A*. Instead, the data consist of observed counts $o_{1},\ldots,o_{K}$ of the number of times allele *A* is observed in *K* different samples of sizes $n_{1},\ldots,n_{K}$ haplotypes, taken at times $t_{1} < \cdots < t_{K}$. For simplicity, we assume that each sampling time *t~k~* is an integer for $k = 1,\ldots,K$. The consecutive observed counts $(o_{k},o_{k + 1},\ldots,o_{k\prime})$ will be denoted by $o_{\lbrack k:k\prime\rbrack}$.

In general, we will denote random variables corresponding to observed quantities using capital letters (e.g., *O~k~*, *C~t~*, and *Y~t~*). The goal is to compute the probability $\mathbb{P}_{\Theta}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ of the observed data, conditional on the model parameters Θ.

Probabilities of Frequency Trajectories
---------------------------------------

Several different evolutionary models can be used to describe stochastic allele frequency changes over time in a population. Discrete changes in allele frequency are often modeled using the Wright--Fisher and Moran processes, whereas continuous changes are often modeled using the diffusion approximation of the Wright--Fisher process ([@msw173-B31]; [@msw173-B13]; [@msw173-B52]) or one of several approximations of the diffusion ([@msw173-B14]; [@msw173-B33]).

Because it is unclear which model provides the most accurate description of biological evolutionary processes, we take the approach in this paper of deriving exact probabilities of allele frequency trajectories under two different evolutionary models: the discrete Wright--Fisher process and the continuous diffusion approximation.

Under the Wright--Fisher model, the probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ of the observed allele counts can be obtained using the recursive formula presented in Procedure 1. Under the diffusion approximation, the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ can be obtained using the recursive formula presented in Procedure 2.

In the "Deterministic allele frequency trajectories under the Wright--Fisher model" and "Deterministic allele frequency trajectories under the diffusion model" sections, we show that if drift is ignored and allele frequencies evolve deterministically, then the probabilities $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ and $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ can be reduced to the simpler approximate probabilities $\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ and $\mathbb{P}_{\Theta,\mathcal{D}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ which ignore the population history and which are computed using Procedures 3 and 4, respectively.

Different estimates of the model parameters Θ can be obtained using each of the different probabilities $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\},\,\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\},\,\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$, and $\mathbb{P}_{\Theta,\mathcal{D}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ by finding the value of Θ that maximizes the given probability of the observed allele counts $o_{\lbrack 1:K\rbrack}$. In our analyses, we estimated the model parameters Θ separately using each of the different probabilities, yielding the estimators ${\widehat{s}}_{\mathcal{W}},\,{\widehat{s}}_{\mathcal{D}},\,{\widehat{s}}_{\mathcal{W}}^{\infty}$, and ${\widehat{s}}_{\mathcal{D}}^{\infty}$. The estimator ${\widehat{s}}_{\mathcal{W}}$ accounts for drift under the discrete Wright--Fisher model, whereas drift in this model is ignored by the estimator ${\widehat{s}}_{\mathcal{W}}^{\infty}$. Similarly, the estimator ${\widehat{s}}_{\mathcal{D}}$ accounts for drift under the diffusion model, whereas drift in this model is ignored by the estimator ${\widehat{s}}_{\mathcal{D}}^{\infty}$.

The degree to which accounting for drift can improve estimates of selection coefficients can be investigated by comparing ${\widehat{s}}_{\mathcal{W}}$ to ${\widehat{s}}_{\mathcal{W}}^{\infty}$ on trajectories simulated under the discrete Wright--Fisher model and by comparing ${\widehat{s}}_{\mathcal{D}}$ to ${\widehat{s}}_{\mathcal{D}}^{\infty}$ on trajectories simulated under the diffusion approximation. **Procedure 1.** Computing $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$1: Define the quantities $\mathbf{d}_{0} = \left( \mathbb{P}\left\{ C_{0} = 0 \right\},\mathbb{P}\left\{ C_{0} = 1 \right\},\ldots,\mathbb{P}\left\{ C_{0} = 2N_{t_{0}} \right\} \right)$ and $\gamma\left( o_{1} \right)$, where $\gamma\left( o_{k} \right) = \left( \gamma_{0}\left( o_{k} \right),\gamma_{1}\left( o_{k} \right),\ldots,\gamma_{2N_{t_{k}}}\left( o_{k} \right) \right)$ with $\gamma_{i}\left( o_{k} \right) = \begin{pmatrix}
n_{k} \\
o_{k} \\
\end{pmatrix}\left( i/2N_{t_{k}} \right)^{o_{k}}\left( 1 - i/2N_{t_{k}} \right)^{n_{k} - o_{k}}$.2: Initialize $\mathbf{v}_{1} = \mathbf{d}_{0}\left\lbrack {\prod\limits_{t = 1}^{t_{1}}\mathbf{T}_{t - 1,t}} \right\rbrack\text{diag}\left\{ \gamma\left( o_{1} \right) \right\}$.3: For $k = 2:K$, compute $$\mathbf{v}_{k} = \mathbf{v}_{k - 1}\left\lbrack {\prod\limits_{t = t_{k - 1} + 1}^{t_{k}}\mathbf{T}_{t - 1,t}} \right\rbrack\text{diag}\left\{ \gamma\left( o_{k} \right) \right\}.$$4: Compute $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = \sum\limits_{i = 0}^{2N_{t_{K}}}v_{K,i}.$\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_**Modifications:** If $\mathbf{d}_{0}$ is unspecified, omit Step 1and set $\mathbf{v}_{1} = \gamma\left( o_{1} \right)/2N_{t_{1}}$ in Step 2. If conditioning on the event *S~K~* that allele *A* is segregating in the final sample, omit Step 4 and instead compute $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}$ using equation(C.1). **Procedure 2.** Computing $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} \right\}$1: For an initial starting frequency *y*~0~ initialize $$\mathbf{b}_{0} = \mathbf{C}_{\ell_{1}}^{- 1}\mathbf{B}_{\ell_{1}}\left( y_{0} \right),$$ where $\mathbf{B}_{\ell}\left( y_{0} \right)$ is the vector of eigenfunctions of the diffusion operator given inequation (A.14) and $\mathbf{C}_{\ell} = \text{diag}\left\{ \left\langle B_{\ell,i},B_{\ell,i} \right\rangle \right\}_{i = 0}^{\infty}$ is given in equation (A.18).2: For $k = 1:K$, compute $$\mathbf{a}_{k}~ = ~\begin{cases}
{\mathbf{b}_{k - 1}\mathbf{E}_{\ell_{k}}\left( t_{k} - t_{k - 1} \right)} & {\,\text{if}\,\ell_{k - 1} = \ell_{k},} \\
{\mathbf{b}_{k - 1}\mathbf{F}\left( t_{k - 1},t_{k};\zeta \right)} & {\,\text{otherwise},} \\
\end{cases}$$ and $$\mathbf{b}_{k} = \mathbf{a}_{k}\mathbf{W}_{\ell_{k}}\mathbf{G}_{\ell_{k}}^{o_{k}}\left( 1 - \mathbf{G}_{\ell_{k}} \right)^{n_{k} - o_{k}}\mathbf{W}_{\ell_{k}}^{- 1},$$ where the matrices $\mathbf{E}_{\ell}\left( t \right),\,\mathbf{F}\left( t_{k - 1},t_{k};\zeta \right),\,\mathbf{W}_{\ell}$, and $\mathbf{G}_{\ell}$ are given byequations (A.17), (B.10), (A.15), and (A.11), respectively, and ζ is the set of Chebyshev nodes in the interval $\left\lbrack 0,1 \right\rbrack$. The matrix inverse $\mathbf{W}_{\ell}^{- 1} = \mathbf{D}_{\ell}\mathbf{W}_{\ell}^{T}\mathbf{C}_{\ell}^{- 1}$ is computed easily using the diagonal matrices $\mathbf{C}_{\ell}$ and $\mathbf{D}_{\ell}$ in equations (A.18) and (A.19).3: Compute $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = \frac{c_{\ell_{K},0}}{B_{\ell_{K},0}\left( 0 \right)}b_{K,0},$$ where $c_{\ell_{K},0} = \left\langle B_{\ell,0},B_{\ell,0} \right\rangle = \left\lbrack C_{\ell_{K}} \right\rbrack_{0,0}$ is the (0, 0) element of matrix $\mathbf{C}_{\ell}$ in equation (A.18) and $B_{\ell_{K},0}\left( 0 \right)$ is the 0^th^ element of the vector $\mathbf{B}_{\ell_{K}}\left( 0 \right)$ in equation (A.14).\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_**Modifications:** If *y*~0~ is unspecified, omit Step 1 and initialize $\mathbf{b}_{1}$ using equation (B.17). Then iterate over $k = 2:K$ in Step 2. If conditioning on the event *S~K~* that allele *A* is segregating in the final sample, omit Step 3 and instead compute $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}$ using equation (D.1). **Procedure 3.** Computing $\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$1: Starting with $y_{0}^{\infty} = y_{0}$, for $t = 0,\ldots,t_{K} - 1$,Compute ${\widetilde{y}}_{t}^{\infty} = u_{aA}^{(t)} + (1 - u_{Aa}^{(t)} - u_{aA}^{(t)})y_{t}^{\infty}$.Compute $$y_{t + 1}^{\infty} = \left\lbrack \frac{{({\widetilde{y}}_{t}^{\infty})}^{2}(1 + s_{t}) + {\widetilde{y}}_{t}^{\infty}(1 - {\widetilde{y}}_{t}^{\infty})(1 + h_{t}s_{t})}{{\overline{w}}_{t}} \right\rbrack,$$ where ${\overline{w}}_{t} = {({\widetilde{y}}_{t}^{\infty})}^{2}(1 + s_{t}) + 2{\widetilde{y}}_{t}^{\infty}(1 - {\widetilde{y}}_{t}^{\infty})(1 + h_{t}s_{t}) + {(1 - {\widetilde{y}}_{t}^{\infty})}^{2}$.2: Compute $$\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = \prod\limits_{k = 1}^{K}{(\begin{matrix}
n_{k} \\
o_{k} \\
\end{matrix})}{(y_{t_{k}}^{\infty})}^{o_{k}}{(1 - y_{t_{k}}^{\infty})}^{n_{k} - o_{k}}.$$\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_**Modifications:** If y~0~ is unspecified, set $\bigtriangleup y = 1/M$ for a large value M and repeat Steps 1 and 2 for the dense uniform grid of M + 1 values $y_{0} \in \lbrack 0, \bigtriangleup y,2 \bigtriangleup y,\ldots,1\rbrack$. Set $\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = \frac{1}{M + 1}\sum\limits_{j = 0}^{M}\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack};y_{0} = j \bigtriangleup y \right\}$. **Procedure 4.** Computing $\mathbb{P}_{\Theta,\mathcal{D}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$1: Fix a large integer *n* and set $\Delta t = 1/n$.2: Starting with $y_{0}^{\infty} = y_{0}$, for $j = 0,\ldots,nt_{K} - 1$, compute $$y_{(j + 1)\Delta t}^{\infty} = {\{{u_{aA}^{(j \bigtriangleup t)} - (u_{aA}^{(j \bigtriangleup t)} + u_{Aa}^{(j \bigtriangleup t)})y_{j\Delta t}^{\infty}}}~~{{+ y_{j\Delta t}^{\infty}(1 - y_{j\Delta t}^{\infty})\lbrack(1 - 2y_{j \bigtriangleup t}^{\infty})h_{j \bigtriangleup t}s_{j \bigtriangleup t} + y_{j\Delta t}^{\infty}s_{j \bigtriangleup t}\rbrack~}\}}\Delta t.$$3: Compute $$\mathbb{P}_{\Theta,\mathcal{D}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = \prod\limits_{k = 1}^{K}{(\begin{matrix}
n_{k} \\
o_{k} \\
\end{matrix})}{(y_{t_{k}}^{\infty})}^{o_{k}}{(1 - y_{t_{k}}^{\infty})}^{n_{k} - o_{k}}.$$\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_**Modifications:** If y~0~ is unspecified, set $\bigtriangleup y = 1/M$ for a large value M and repeat Steps 1 and 2 for the dense uniform grid of M + 1 values $y_{0} \in \lbrack 0, \bigtriangleup y,2 \bigtriangleup y,\ldots,1\rbrack$. Set $\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = \frac{1}{M + 1}\sum\limits_{j = 0}^{M}\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack};y_{0} = j \bigtriangleup y \right\}$.

Overview of the Experimental Design
-----------------------------------

We simulated allele frequency trajectories under a variety of selection strengths and piecewise constant population histories reflecting demographic patterns such as exponential growth, bottlenecks, rapid population size oscillations, and constant histories. We then compared the demography-aware estimates ${\widehat{s}}_{\mathcal{W}}$ and ${\widehat{s}}_{\mathcal{D}}$ with the estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$ and ${\widehat{s}}_{\mathcal{D}}^{\infty}$ that ignore drift to study the degree to which accounting for population size can improve the accuracy of inferences.

Expected Allele Frequency Trajectories
--------------------------------------

Before comparing the accuracy of the different estimators, we first explored the degree to which trajectories that ignore drift differ from trajectories that account for drift resulting from finite population sizes. [Figure 2](#msw173-F2){ref-type="fig"} shows the expected frequency of allele *A* in a discrete Wright--Fisher population of constant size for several different initial allele frequencies, selection coefficients, and effective population sizes. [Figure 2](#msw173-F2){ref-type="fig"} illustrates that, for any starting frequency and selection coefficient, the mean allele frequency trajectory approaches the mean trajectory in a population without drift (i.e., in a population of infinite size), as the true population size increases. Moreover, if the initial frequency is sufficiently high, the expected trajectory is close to its deterministic limit even when the population size is small and drift is high. Fig. 2.Expected Wright--Fisher trajectories of allele *A* for different initial starting counts *c*~0~, selection coefficients *s*, and effective population sizes *N*. Columns correspond to different initial starting frequencies $c_{0}/2N$ with $c_{0}/2N = 0.01,0.1,$ and 0.25. The dominance parameter is set to *h* = ½ in all panels. Because the effects of mutation are negligible during the time periods we consider, we set $u_{Aa} = u_{aA} = 0$.

From [figure 2](#msw173-F2){ref-type="fig"}, it can be seen that an effective population size of several thousand individuals is often sufficiently large to produce deterministic behavior, even when the selection coefficient and initial allele frequency are small. Thus, selection coefficient inference methods that ignore drift are likely to be accurate for a broad range of population sizes and selection coefficients. As we will see, methods that ignore drift can be almost as accurate as methods that account for drift, even within the small-parameter-value regime.

Inference Accuracy for Populations of Constant Size
---------------------------------------------------

To explore how accounting for drift affects inference accuracy, we first considered the accuracy of inferring selection coefficients in a population of constant finite size. [Figure 3](#msw173-F3){ref-type="fig"} shows the maximum likelihood estimate (MLE) of the selection coefficient for three different effective population sizes (*N* = 100, 500, 1000), three selection coefficients $(s = 0.01,0.05,0.1)$, and two initial allele frequencies $(y_{0} = 0.01,0.1)$ for *h* = ½. In each of panels A-R, the violin plots summarize the maximum likelihood estimates for 100 different simulation replicates in which an allele frequency trajectory was simulated for 500 generations with samples of size *n* = 50 taken at generations $t =$50, 100, 150, 200, 250, 300, 350, 400, 450, and 500. Fig. 3.Maximum likelihood estimates of the selection coefficient *s* in populations of constant size. For each of three different selection coefficients ($s = 0.01,0.05,0.1$) and effective population sizes (*N* = 100, 500, 1000), 100 allele frequency trajectories were simulated for 500 generations under the either the Wright--Fisher or diffusion models. Samples of 50 alleles were taken at times 50, 100, 150, 200, 250, 300, 350, 400, 450, and 500 generations. Bars in panels A-R indicate medians. Dashed lines indicate true selection coefficients. The maximum width of each violin plot is scaled to the same value for all estimators. Scatter plots S, T, U, V, W, and X compare the estimates ${\widehat{s}}_{\mathcal{W}}$ with the estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$ using the same data shown in the top modes of the distributions in Panels A, B, C, J, K, and L for the case *s* = 0.01. Diagonal lines in these scatter plots indicate the line ${\widehat{s}}_{\mathcal{W}}^{\infty} = {\widehat{s}}_{\mathcal{W}}$. Bimodal violin plots are due to the fact that allele frequency trajectories typically fall into one of two categories: trajectories in which allele *A* is lost quickly, resulting in a strong negative estimate of the selection coefficient, and trajectories in which allele *A* remains segregating long enough to allow a more accurate estimate of the selection coefficient.

For the discrete Wright--Fisher model, allele frequency trajectories were simulated by sampling the allele frequency in each generation from the vector of transition probabilities, conditional on the frequency in the previous generation using Procedure 5. Under the diffusion model, trajectories were sampled using the approach in Procedure 6. Maximum likelihood estimates were obtained for the Wright--Fisher trajectories using the grid search described in Procedure 7 over the likelihoods computed using Procedures 1 and 3, and maximum likelihood estimates for the diffusion trajectories were obtained using the same grid search approach over the likelihoods computed using Procedures 2 and 4. In each panel in [figure 3](#msw173-F3){ref-type="fig"}, the estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$ and ${\widehat{s}}_{\mathcal{W}}$ were computed for the same set of 100 allele frequency trajectories simulated under the discrete Wright--Fisher model and the estimates ${\widehat{s}}_{\mathcal{D}}^{\infty}$ and ${\widehat{s}}_{\mathcal{D}}$ were computed for the same set of 100 allele frequency trajectories simulated under the diffusion model.

By comparing the estimates computed accounting for drift with the estimates obtained ignoring drift, it can be seen that all methods have similar accuracies. All methods perform well when the population size, selection coefficient, and initial frequency are sufficiently large (e.g., [fig. 3I](#msw173-F3){ref-type="fig"} for the case $y_{0} = 0.01$ and Panels 3K through 3R for the case $y_{0} = 0.1$), and all methods have reduced accuracy, otherwise. To put this another way: the allele frequency trajectories for which selection coefficients are inferred accurately by demography-aware methods correspond to those for which the deterministic estimates are also accurate. Thus, methods that ignore or account for drift are likely to produce estimates with similar accuracy.

Moreover, it can be seen from the scatter plots (Panels 3S--3X), which compare the estimators ${\widehat{s}}_{\mathcal{W}}$ and ${\widehat{s}}_{\mathcal{W}}^{\infty}$, that the point estimates themselves are very similar for both the demography-aware and deterministic methods. Although this similarity may be expected given that the deterministic methods differ from the demography-aware methods only in that they ignore the additional variability in the allele frequencies arising from genetic drift, it is surprising that the point estimates are so similar, as the overall expected allele frequency trajectory in the deterministic case can differ considerably from the expected trajectory accounting for drift ([fig. 2A*--C*](#msw173-F2){ref-type="fig"}).

As the magnitude of the selection coefficient decreases, the point estimates of the deterministic and demography-aware estimators remain well correlated, although the accuracy of all methods decreases. This can be seen in [figure 3S--*X*](#msw173-F3){ref-type="fig"} and in [supplementary fig. S1, Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msw173/-/DC1) online in which the estimates by the different methods remain correlated, but become more variable as the selection coefficient *s* decreases and an increasingly large number of trajectories drift out of the population quickly, leading to strong negative estimates of selection coefficients by both methods.

Inference Accuracy in Populations of Piecewise Constant Size
------------------------------------------------------------

We next explored the degree to which accounting for more complicated population histories can improve maximum likelihood estimates, focusing on three scenarios, a population with a bottleneck, a population undergoing exponential growth, and a population undergoing rapid oscillations in size. Under each scenario, we simulated 100 allele frequency trajectories for an allele with selection coefficient *s* = 0.05, dominance parameter *h* = ½, and initial frequency $y_{0} = 0.1$ under the Wright--Fisher model and separately under the diffusion model. The parameter values in these simulations were chosen so that drift would be strong enough to affect allele frequency trajectories, but not strong enough to result in poor estimates of selection coefficients by the full-likelihood methods.

In addition to comparing estimates made by the deterministic estimators, ${\widehat{s}}_{\mathcal{W}}^{\infty}$ and ${\widehat{s}}_{\mathcal{D}}^{\infty}$, with those of the exact estimators, ${\widehat{s}}_{\mathcal{W}}$ and ${\widehat{s}}_{\mathcal{D}}$, that account for the true time-varying population history, we investigated the effect on accuracy of using crude, yet reasonable estimates of the population history. In particular, we also inferred selection coefficients using likelihoods computed using variants of Procedures 1 and 2 in which the population was assumed to consist of a single epoch of constant size equal to the Watterson estimate ([@msw173-B54], Eqn. 1.4a; [@msw173-B27], p.62). The Watterson estimate was obtained by computing the expected site frequency spectrum (SFS) for the multi-epoch model for a sample size of 20 alleles using the method of [@msw173-B30] and then inferring the effective size of a single epoch using Watterson's estimator (Computing the Watterson estimator of *N*~e~ from the expected SFS of a piecewise constant population). The discrete Wright--Fisher and diffusion estimators based on the Watterson estimate of effective size are denoted by ${\widehat{s}}_{\mathcal{W}}^{N_{e}}$ and ${\widehat{s}}_{\mathcal{D}}^{N_{e}}$, respectively.

### The Case of a Bottleneck

To model populations with bottlenecks, we considered populations composed of three epochs, each of length 100 generations, with sizes *N*~1~, *N*~2~, and *N*~3~ satisfying $N_{1} = N_{3} = 5N_{2}$. Samples of size 50 were taken at times $50,\, 100,\, 150,\, 200,\, 250$, and 300. [Figure 4A and B](#msw173-F4){ref-type="fig"} shows the results for two different populations; in the population in [figure 4A](#msw173-F4){ref-type="fig"}, we set $N_{1} = 500$ and in the population in [figure 4B](#msw173-F4){ref-type="fig"} we set $N_{1} = 2500$. Fig. 4.Maximum likelihood estimates of the selection coefficient *s* in populations with a bottleneck, exponential growth, or rapidly oscillating population size. In each panel, the trajectory of an allele with selection coefficient *s* = 0.05, dominance parameter *h* = ½, and starting frequency $y_{0} = 0.1$ was simulated 100 times under the Wright--Fisher and diffusion models. Bars indicate medians. Dashed lines indicate the true selection coefficient. The maximum width of each violin plot is scaled to the same value for all estimators.

From [figure 4A and *B*](#msw173-F4){ref-type="fig"}, it can be seen that all methods performed similarly. However, the deterministic estimators had significantly lower bias than the full-likelihood estimators computed using the mis-specified population history for a bottleneck of size $N_{2} = 100$ with $N_{1} = N_{3} = 500$ ([fig. 4A](#msw173-F4){ref-type="fig"}). In the case of the bottleneck in [figure 4A](#msw173-F4){ref-type="fig"}, the means of the deterministic estimators ${\widehat{s}}_{\mathcal{W}}^{\infty}$ and ${\widehat{s}}_{\mathcal{D}}^{\infty}$ were not significantly different from the true selection coefficient *s* (*P* = 0.59 and *P* = 0.93 for *t*-tests of the null hypotheses $\text{mean}({\widehat{s}}_{\mathcal{W}}^{\infty}) = s$ and $\text{mean}({\widehat{s}}_{\mathcal{D}}^{\infty}) = s$ versus the alternative hypotheses $\text{mean}({\widehat{s}}_{\mathcal{W}}^{\infty}) \neq s$ and $\text{mean}({\widehat{s}}_{\mathcal{D}}^{\infty}) \neq s$). In comparison, the means of the full likelihood estimators ${\widehat{s}}_{\mathcal{W}}^{N_{e}}$ and ${\widehat{s}}_{\mathcal{D}}^{N_{e}}$ with mis-specified histories were significantly different from *s* (*P* = 0.03 and *P* \< 0.01, respectively). A similar trend was observed for the bottleneck history with larger sizes shown in [figure 4B](#msw173-F4){ref-type="fig"} (*P* = 0.64 and *P* = 0.12 for *t*-tests of the null hypotheses $\text{mean}({\widehat{s}}_{\mathcal{W}}^{\infty}) = s$ and $\text{mean}({\widehat{s}}_{\mathcal{D}}^{\infty}) = s$, compared with *P* = 0.13 and *P* = 0.01 for *t*-tests of the null hypotheses $\text{mean}({\widehat{s}}_{\mathcal{W}}^{N_{e}}) = s$ and $\text{mean}({\widehat{s}}_{\mathcal{D}}^{N_{e}}) = s$). Note that, despite the tight bottleneck in [figure 4A](#msw173-F4){ref-type="fig"}, inferences were still relatively accurate due to the larger sizes of epochs 1 and 3.

### The Case of Exponential Growth

To model exponential growth, we considered populations composed of five epochs, each of length 100 generations, with effective population sizes chosen to represent 5-fold exponential growth across all five epochs. Specifically, the size in epoch ℓ was set to $N_{\ell} = N_{0}e^{\eta\tau_{\ell}}$, where we chose $N_{0} = 100$ and the growth constant *η* was chosen such that $e^{\eta\tau_{5}} = 5$. Samples of size 50 were taken in generations 100, 200, 300, 400, and 500. From the results in [figure 4C and D](#msw173-F4){ref-type="fig"}, it can be seen that all methods performed with similar accuracy in the growth scenario.

### The Case of Rapidly Oscillating Population Size

[Figure 4E and F](#msw173-F4){ref-type="fig"} shows inferences of the selection coefficient in a population with rapidly oscillating size. Such demographic histories, which are often seen in insect populations like *Drosophila*, have moderate arithmetic mean sizes, but small harmonic mean sizes and experience episodes of extreme drift.

In the simulations shown in [figure 4E](#msw173-F4){ref-type="fig"}, the population size oscillates rapidly between 10 and 100 diploids every five generations. In the simulations shown in [figure 4F](#msw173-F4){ref-type="fig"}, the population size oscillates between 100 and 500 diploids every five generations. From [figure 4](#msw173-F4){ref-type="fig"}, it can be seen that the methods that ignore drift have similar accuracy to the methods that account for drift. However, the methods that account for drift are slightly less biased when the population size oscillates between very small values ([fig. 4E](#msw173-F4){ref-type="fig"}) (*P* = 0.017 and *P* = 0.012, respectively, for the one-tailed *t*-tests of the null hypotheses ${\widehat{s}}_{\mathcal{W}}^{\infty} \geq {\widehat{s}}_{\mathcal{W}}$ and ${\widehat{s}}_{\mathcal{D}}^{\infty} \geq {\widehat{s}}_{\mathcal{D}}$ versus the alternative hypotheses ${\widehat{s}}_{\mathcal{W}}^{\infty} < {\widehat{s}}_{\mathcal{W}}$ and ${\widehat{s}}_{\mathcal{D}}^{\infty} < {\widehat{s}}_{\mathcal{D}}$).

Conditioning on Segregation in the Final Sample
-----------------------------------------------

It is sometimes of interest to infer the selection coefficient of an allele, conditional on the event that the allele is segregating in the most recent sample. Such conditional inferences are useful if alleles are ascertained in present-day samples and their historical trajectories are subsequently investigated.

Conditioning on segregation in the final sample is also useful for estimating weak positive selection coefficients when initial allele frequencies are low. This is because a large fraction of weakly selected alleles with low initial frequencies will drift out of the population quickly resulting in large negative estimates of their selection coefficients. However, more accurate estimates can be obtained for the subset of alleles that are not lost quickly, which can be seen, for example, in [figure 3B--*I*](#msw173-F3){ref-type="fig"} through I in which the part of the density corresponding to alleles that are not lost quickly from the population is localized around the true selection coefficient.

Considering only alleles that are segregating in the final sample can lead to biased estimates of selection coefficients if likelihood methods do not properly condition on segregation. For example, weakly selected alleles typically drift out of small populations quickly. Thus, weakly selected alleles that escape loss and ultimately fix generally exhibit faster-than-expected increases in frequency that are similar to the unconditional trajectories of alleles under stronger selection. Thus, if a likelihood method does not properly account for conditioning, weakly selected alleles that are segregating in the final sample will have inflated inferred selection coefficients.

Estimators that ignore drift cannot be modified to condition on the event of segregation in the final sample because they implicitly assume that alleles follow fixed trajectories whose long-term behavior in the absence of mutation is entirely determined by the selection coefficient: fixation for positively selected alleles and loss for negatively selected alleles. Thus, estimators that ignore drift are expected to produce biased estimates of selection coefficients when applied to conditioned trajectories.

In contrast, the allele frequency trajectories in likelihood methods that account for the population size are modeled stochastically, allowing likelihoods to be modified to condition on segregation in the final sample. It is expected that methods that account for the true population size can be modified to produce accurate estimates of selection coefficients, whereas methods that ignore drift will necessarily produce biased estimates.

### Simulations Conditioning on Segregation

To investigate the degree to which accounting for drift can improve estimates of selection coefficients when allele frequency trajectories are conditioned on segregation in the final sample, we modified the discrete Wright--Fisher probability in Procedure 1 to compute the likelihood conditional on segregation in the final sample using results derived in the "Conditional Probabilities" section. Under a grid search, this modified likelihood yields the conditional maximum likelihood estimator ${\widehat{s}}_{\mathcal{W}|S_{K}}$. We compared the estimates computed using the exact conditional estimator ${\widehat{s}}_{\mathcal{W}|S_{K}}$ with estimates computed using the approximate estimator ${\widehat{s}}_{\mathcal{W}}^{\infty}$ that ignores drift and cannot be modified to account for conditional allele frequency trajectories.

The effect of failing to account for conditioning is evident in the violin plots in [figure 5A--*I*](#msw173-F5){ref-type="fig"} corresponding to the unconditional approximate maximum likelihood estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$. As expected, when the true selection coefficient is small ($s \leq 0.01$), the estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$ are biased upward. Conversely, when the selection coefficient is larger ($s \geq 0.05$), the approximate estimator ${\widehat{s}}_{\mathcal{W}}^{\infty}$ produces negatively biased estimates because alleles under strong positive selection that remain segregating in the final sample show slower-than-expected increases in frequency. In contrast to the estimator ${\widehat{s}}_{\mathcal{W}}^{\infty}$, the bias is negligible in the estimator ${\widehat{s}}_{\mathcal{W}|S_{K}}$, which accounts for drift and properly conditions on segregation in the final sample. Fig. 5.Estimates of selection coefficients, conditional on segregation. Each violin plot was computed using 100 frequency trajectories sampled over 500 generations for an allele with initial frequency $y_{0} = 0.01$. As in [figure 3](#msw173-F3){ref-type="fig"}, samples of size *n* = 50 were taken in generations 50, 100, 150, 200, 250, 300, 350, 400, 450, and 500. In Panels A--I, trajectories were sampled conditional on the event that the selected allele was segregating in the final sample. In Panels J--R, trajectories were sampled conditional on the event that the selected allele was either segregating or fixed in the final sample. Red bars indicate medians. Dashed lines indicate the true selection coefficients. The maximum width of each violin plot is scaled to the same value for both estimators.

The results shown in [figure 5A--*I*](#msw173-F5){ref-type="fig"} suggest that methods that account for drift are capable of significantly improving the accuracy of estimates of selection coefficients when allele frequency trajectories are conditioned on segregation. The differences in accuracy between methods that ignore or account for drift are visible for a range of selection coefficients and population sizes. However, the differences in accuracy between the methods diminish as the population size becomes larger.

### Simulations Conditioning on Segregation or Fixation

The magnitude of the bias in the estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$ is due in part to the event on which trajectories are conditioned. In cases involving positive selection in populations of moderate or large size, most alleles will be fixed in the final sample (e.g., $> 80\%$ fixation within 10 generations when *s* = 0.1, *h* = ½, $y_{0} = 0.01$, and *N* = 1000). Thus, it may sometimes be more natural to condition on the event *F~K~* that a selected allele is found (segregating or fixed) in the final sample. Under this conditioning scheme, the approximate estimator ${\widehat{s}}_{\mathcal{W}}^{\infty}$ will not generally produce negatively biased estimates of selection coefficients because allele frequency trajectories will not be constrained to those which exhibit slower-than-expected increases in allele frequency.

In light of these considerations, we repeated the analysis shown in [figure 5A--*I*](#msw173-F5){ref-type="fig"}, simulating allele frequency trajectories conditional on the event that the allele was segregating or fixed in the final sample. To compare the estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$ with maximum likelihood estimates that fully account for drift and the proper conditioning, we also modified the probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ computed in Procedure 1 to condition on the event *F~K~* of segregation or fixation in the final sample, yielding the conditional probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| F_{K} \right\}$ ([eq. 21](#E29){ref-type="disp-formula"}) with the associated estimator ${\widehat{s}}_{\mathcal{W}|F_{K}}$.

By comparing [figure 5J--*R*](#msw173-F5){ref-type="fig"} with [figure 5A--*I*](#msw173-F5){ref-type="fig"}, it can be seen that the estimator ${\widehat{s}}_{\mathcal{W}}^{\infty}$ has considerably less bias when conditioning on the event *F~K~* than when conditioning on *S~K~*. Although the bias is still high when the population size is small ($N \approx 100$), it decreases quickly as the population size increases and becomes comparable to the bias in the properly conditioned, demography-aware estimator ${\widehat{s}}_{\mathcal{W}|F_{K}}$ when the population size is greater than approximately *N* = 500 diploids. In contrast to [figure 5E--*I*](#msw173-F5){ref-type="fig"}, the bias in ${\widehat{s}}_{\mathcal{W}}^{\infty}$ observed in [figure 5M--*R*](#msw173-F5){ref-type="fig"} is positive because the trajectories on which these estimates are based exclude only those in which the allele is lost; thus, they exhibit faster-than-expected growth on average. The results in [figure 5J--*R*](#msw173-F5){ref-type="fig"} suggest that under certain conditioning schemes, methods that ignore drift can produce similar estimates to methods that account for drift.

The Effect of Sample Size on Accuracy
-------------------------------------

When the sample size is small, the variance in estimates arising from sampling noise will tend to obscure small differences between estimators that ignore or account for population size. Thus, when comparing methods, it is important to sample a sufficiently large number of alleles to ensure that the differences between the methods due to ignoring or accounting for drift are visible.

To evaluate the effects of sample size on inference accuracy, we inferred the selection coefficient for a range of sample sizes for several different combinations of the population size and selection coefficient. [Figure 6](#msw173-F6){ref-type="fig"} shows a plot of the variance in selection coefficients inferred using Procedures 1 and 3 for sample sizes ranging from *n* = 2 to *n* = 50. For each combination of *N~e~*, *s*, and *n*, the trajectories of 100 alleles were simulated under the Wright--Fisher process with an initial allele frequency of $y_{0} = 0.1$. Samples were taken in generations 50 and 100. Fig. 6.The effect of sample size on inference accuracy. The variance of the estimates produced by the methods in Procedures 1 and 3 are shown for a range of sample sizes.

The plots in [figure 6](#msw173-F6){ref-type="fig"} suggest that variability due to small sample sizes has a strong effect on the variability in estimates only for sample sizes smaller than 10 alleles. Thus, in all of our simulations, we have used a sample size of *n* = 50 alleles so that differences between estimators are not likely to be obscured by the variance in estimates due to small sample sizes.

Unspecified Initial Allele Frequencies
--------------------------------------

In our simulations, we have assumed that the initial frequency *y*~0~ of allele *A* at time *t* = 0 is known. Knowledge of the initial allele frequency provides approximately the same information as an informative sample taken at time *t* = 0. Although this information can improve the estimate of the selection coefficient if few samples are taken, it contributes relatively little information to the estimate of *s* when the number of samples is moderate. However, because a lack of knowledge of the initial frequency could potentially lead to greater errors in the deterministic estimators, we also compared versions of the estimators that assume a uniform prior distribution on the allele frequency at the time of the first sampling event, and which do not incorporate information about the initial allele frequency (Computing $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ When $\mathbf{d}_{0}$ Is Unspecified, Computing $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ When *y*~0~ Is Unspecified, and Deterministic Estimators When the Initial Allele Frequency Is Unknown).

[Figure 7](#msw173-F7){ref-type="fig"} shows a comparison between the deterministic and exact Wright--Fisher estimates when the initial frequency is drawn uniformly from the interval $y_{0} \in (0,0.5\rbrack$ and the selection coefficient is inferred using Procedures 1 and 2 assuming no knowledge of the initial allele frequency. From [figure 7](#msw173-F7){ref-type="fig"}, it can be seen that the deterministic and exact estimators produce similar results, even in the absence of knowledge of the initial frequency. Moreover, although the estimates produced by the two methods decrease in accuracy as the number of sampling events decreases, the accuracies of the two kinds of methods remain similar. In particular, From [figure 7D--*F*](#msw173-F7){ref-type="fig"}, it can be seen that both the deterministic and exact methods are relatively accurate when allele frequency trajectories do not drift out or fix immediately due to random fluctuations. Fig. 7.Scatter plots comparing the estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$ and ${\widehat{s}}_{\mathcal{W}}$ when the initial allele frequency is unspecified during inference. One hundred allele frequency trajectories were simulated for 500 generations with *s* = 0.05, *N* = 100, and *y*~0~ sampled uniformly from the interval $(0,0.5\rbrack$ as described in the "Simulations" section. Estimates ${\widehat{s}}_{\mathcal{W}}^{\infty}$ and ${\widehat{s}}_{\mathcal{W}}$ were computed using Procedures 1 and 2 using the modified procedure taken when $\mathbf{d}_{0}$ and *y*~0~ are unspecified. (A) Estimates computed using all ten samples of size *n* = 50 taken in generations $50,100,150,200,250,300,350,400,450$, and 500. (B) Estimates computed using the five samples taken in generations $50,150,250,350$, and 450. (C) Estimates computed using the two samples taken in generations 50 and 150. (D), (E), and (F) are zoomed-in versions of the plots in (A), (B), and (C). All panels correspond to the same 100 simulated trajectories. Dashed lines correspond to the true selection coefficient. Bars show means over all 100 estimates. The solid line is ${\widehat{s}}_{\mathcal{W}}^{\infty} = {\widehat{s}}_{\mathcal{W}}$.

Violation of Model Assumptions
------------------------------

It is possible that violations of model assumptions in real data could increase the differences in accuracy between the deterministic and demography-aware estimators. To investigate this possibility we compared the performance of the deterministic and demography-aware discrete Wright--Fisher estimators using experimentally sampled allele frequency time series data that are thought to violate several model assumptions. In particular, we considered time series data for an allele at the *medionigra* locus in the species *Panaxia dominula*, which confers a darkened wing phenotype ([@msw173-B7]). Although these samples are thought to represent a single population that is isolated from migration events, there is evidence for temporal fluctuations in the selection coefficient, frequency dependent selection, and assortative mating. These data are also useful for our analyses because they include estimates of population sizes obtained using mark, release, and recapture. The data set, which spans over 40 years, provided important evidence for natural selection resulting from environmental pressures.

We computed the log likelihood as a function of the selection coefficient *s* for three different values of the dominance parameter ($h = 0,0.5$, and 1). [Figure 8](#msw173-F8){ref-type="fig"} shows the log likelihood for different values of *h*, along with asymptotic normal approximations of 95% confidence intervals. [Figure 8](#msw173-F8){ref-type="fig"} shows that the log likelihood surfaces computed using Procedures 1 and 3 are qualitatively similar, yielding similar point estimates for *s*. However, the confidence intervals for the deterministic estimator are considerably smaller than those of the demography-aware estimator. This result is expected, given that the deterministic estimator ignores the largely symmetrical variability arising from drift and considers only variability in the sampling frequency. Fig. 8.The log likelihood as a function of *s* computed using Procedures 1 and 3 for three different values of the dominance parameter *h* ($h = 0,0.5,1$). Solid vertical lines show the maximum likelihood estimates, ${\widehat{s}}_{\mathcal{W}}$ in (A) and ${\widehat{s}}_{\mathcal{W}}^{\infty}$ in (B), for the curve of the corresponding color. Dashed vertical lines indicate asymptotic normal approximations of the 95% confidence intervals for *s*, for the curve of the corresponding color. (A) The demography-aware likelihoods (Procedure 1). (B) The deterministic likelihoods (Procedure 3).

Note that the point estimates presented in [figure 8](#msw173-F8){ref-type="fig"} are similar to those estimated previously in other studies. Using a model of additive selection (*h* = 0.5) [@msw173-B7] estimated $s \approx - 0.16$ whereas [@msw173-B38] inferred $s \approx - 0.12$. These estimates are close to the maximum likelihood estimates that we inferred under the same model of additive selection (${\widehat{s}}_{\mathcal{W}} = - 0.1$ and ${\widehat{s}}_{\mathcal{W}}^{\infty} = - 0.09$) (note that the estimates in [@msw173-B7] and [@msw173-B38] are reported using a different parameterization of the dominance model than the one we have used in this paper; thus, we have scaled the selection coefficients reported in these papers so that they are comparable with the ones reported here).

It is of interest to note that [@msw173-B38] found evidence that the *medionigra* allele is recessive, finding that the likelihood was maximized for a recessive model with a strong negative selection coefficient around $s \approx - 1$, although they note that such a strong negative selection coefficient violates the assumptions of the Gaussian model under which their likelihoods were derived. In accordance with the result of [@msw173-B38], we find that both the deterministic and exact Wright--Fisher likelihoods we considered were also maximized for a recessive pattern of dominance (*h* = 0) with a large negative selection coefficient around $s \approx - 1$.

Under a recessive model of dominance, the strongly negative selection coefficient that was inferred by the estimators in this study and by that of [@msw173-B38] is reasonable given that the allele frequency decreases rapidly, which would be unlikely under weak selection where alleles must combine in homozygotes for selection to act. Unlike the analyses of [@msw173-B38], a selection coefficient of $s \approx - 1$ does not violate our model assumptions. Thus, our analysis provides further support that the *medionigra* variant is largely recessive, although there is evidence that the dominance of the variant can change over time ([@msw173-B7]).

Computational Efficiency of Methods
-----------------------------------

As we have noted, methods that assume that allele frequency trajectories are deterministic can be considerably faster than methods that account for population size histories. [Table 1](#msw173-T1){ref-type="table"} shows the average runtimes of the estimators ${\widehat{s}}_{\mathcal{W}}^{\infty},\,{\widehat{s}}_{\mathcal{W}},\,{\widehat{s}}_{\mathcal{D}}$, and ${\widehat{s}}_{\mathcal{D}}^{\infty}$ for the computations used to produce [figure 3A--*I*](#msw173-F3){ref-type="fig"}. Table 1Mean Runtimes of the Methods in [figure 3A--*I*](#msw173-F3){ref-type="fig"} (seconds).*N*~e~s${\widehat{s}}_{\mathcal{W}}^{\infty}$${\widehat{s}}_{\mathcal{W}}$${\widehat{s}}_{\mathcal{D}}^{\infty}$${\widehat{s}}_{\mathcal{D}}$1000.010.012.304.74197.250.050.012.364.23204.660.10.012.293.98217.345000.010.02134.074.41185.180.050.01132.454.41496.830.10.02126.354.46531.1510000.010.02175.274.64196.900.050.02191.534.78815.270.10.02199.324.671950.59

From the table, it can be seen that the runtimes are considerably faster for the estimators based on deterministic trajectories (${\widehat{s}}_{\mathcal{W}}^{\infty}$ and ${\widehat{s}}_{\mathcal{D}}^{\infty}$). Moreover, the runtimes for ${\widehat{s}}_{\mathcal{W}}^{\infty}$ and ${\widehat{s}}_{\mathcal{D}}^{\infty}$ do not depend on the population size or selection coefficient. In comparison, the runtimes for the estimators ${\widehat{s}}_{\mathcal{W}}$ and ${\widehat{s}}_{\mathcal{D}}$ increase with increasing *N~e~* and *s* because these methods depend on eigenvalue decompositions or sparse matrix products, which require larger matrices or greater precision when *N~e~* or *s* is large. Note that although the discrete Wright--Fisher estimator is considerably faster than the diffusion estimator for the scenarios we considered, the diffusion estimator can still be more efficient when samples are widely separated in time and the repeated matrix--vector products required by the discrete Wright--Fisher method become cumbersome.

Discussion
==========

The results of our analyses suggest that accurate estimates of selection coefficients from allele frequency time series data can often be obtained by assuming that alleles evolve without drift in a population of infinite size. In the majority of our simulations, the estimates obtained using deterministic approximations were nearly as accurate as estimates obtained by explicitly modeling the true population history and they were sometimes more accurate than estimates obtained using crude but reasonable estimates of the population history. The latter result indicates that an overestimate of the population size, for example due to underestimation of the mutation rate or other factors, may have only small adverse effects on the accuracy of exact likelihood estimates of the selection coefficient from time series data, a result that is of interest because population size histories can be difficult to infer accurately.

Surprisingly, estimates made under the deterministic approximation were generally as accurate as estimates that accounted for drift, due to the fact that the exact maximum likelihood methods had low accuracy when drift was strong. Accounting for the true population history only resulted in significantly improved estimates of selection coefficients when conditioning on the event that the target allele was segregating in the final sample. Methods that modeled the true population history were more accurate in this case because they could be modified to model conditional trajectories, whereas methods that assumed infinite population sizes could not. These results suggest that methods that account for drift are likely to be preferable under circumstances in which conditioning on segregation is desirable. However, it is important to note that deterministic methods can perform well when population sizes are moderately large if allele frequencies are conditioned on a slightly different event: the event that an allele is found (segregating or fixed) in the final sample.

It is important to note that we have focused on estimates of the selection coefficient rather than the dominance parameter *h*, which is a more difficult task. However, in our analyses of the *medionigra* variant in the *Panaxia dominula* moth, we found that both the deterministic and exact likelihoods were qualitatively similar in both *h* and *s* and were both maximized for a recessive model of dominance. A comprehensive exploration of the effects of deterministic assumptions on the inference of *h* is deferred to future analyses.

It is also important to note that our analyses do not imply that the deterministic and exact estimators will produce exactly the same point estimate for a given allele frequency trajectory. Although the scatter plots in [figures 3, 7](#msw173-F3 msw173-F7){ref-type="fig"}, and [supplementary figure S1, Supplementary Material](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msw173/-/DC1) online suggest that both deterministic and exact estimators often agree closely on the magnitude of the selective strength, the estimates for any given trajectory can differ by 100% or more. This result is consistent with the observation by [@msw173-B45] that specifying different demographic histories led to different estimates of dominance and selection parameters for alleles affecting coat coloration in horses. Although specific point estimates can differ, our results demonstrate that deterministic and exact estimators have similar accuracies for estimating the strength of selection, even when the population size is small.

The idea that ignoring drift can lead to accurate estimates of selection coefficients is not new. In fact, inference methods based on deterministic allele frequency trajectories capitalize on exactly this idea. However, our comparison with estimators based on exact likelihoods makes it possible to characterize the relative loss in accuracy that is incurred when drift is ignored, as well as the demographic, evolutionary, and sampling scenarios under which accounting for drift is likely to be important.

The comparatively accurate estimates achieved by methods that assume deterministic allele frequency trajectories are encouraging for three primary reasons. First, a large number of studies have relied on the assumption that alleles evolve deterministically in order to infer selection coefficients from biological time series data. Our results suggest that these estimates are likely to be nearly as accurate as those obtained using the exact likelihood accounting for drift. Second, estimators based on deterministic trajectories can be considerably faster than estimators that account for drift, making them useful for inferring selection coefficients at large numbers of loci. Third, it may be easier to obtain analytical results under the assumption that allele frequencies change deterministically, simplifying the development of inference methods for inferring selection coefficients under more complicated scenarios; for example, inferring coefficients at linked loci ([@msw173-B28]). Finally, the ability to ignore the population size is useful in situations in which the true population history is unknown or difficult to infer.

In addition to characterizing the degree to which accounting for drift can improve estimates of selection coefficients, our results shed light on the accuracy of exact maximum likelihood methods for inferring selection coefficients from allele frequency trajectories. In accordance with predictions about the relative strengths of genetic drift and selection ([@msw173-B21], section 3.7) and experimental work ([@msw173-B20]), our findings suggest that very small selection coefficients ($s \leq 0.01$) are difficult to infer if the initial allele frequency and population size are not large. Moreover, even if the population size is large, the accurate inference of a small selection coefficient may require samples taken over hundreds of generations, during which time the selection coefficient could change considerably ([@msw173-B15]; [@msw173-B46]).

Despite the difficulties of inferring weak selection coefficients when the population size is small, coefficients of one percent or lower can be inferred accurately if the initial allele frequency is sufficiently high. It is important to note that the selection coefficient need not be high at the time of the very first sampling event, as long as the allele has reached a sufficiently high frequency at one of the intermediate sampling events, leading to quasi-deterministic behavior between some sampling time points that can be exploited by the maximum likelihood estimator.

Although we have only considered positively selected alleles in our simulation analyses, our results apply equally well to negatively selected alleles, as it is arbitrary whether we choose to track the trajectory of the allele with higher or lower fitness. We have also focused on low initial allele frequencies $(y_{0} \leq 0.1)$ for selected alleles; however, it is clear from [figure 2](#msw173-F2){ref-type="fig"} that allele frequency trajectories become increasingly deterministic as the initial allele frequency increases. Thus, the accuracy of a method that assumes a deterministic trajectory will become more similar to that of a method that accounts for drift as the initial allele frequency increases. Conversely, for negatively selected alleles, the accuracy of the deterministic method will approach that of the exact likelihood as the initial allele frequency decreases. Thus, our analyses provide a characterization of inference accuracy for both positively and negatively selected alleles for a broad range of starting frequencies.

At first glance, our finding that the population size does not strongly influence estimates of selection coefficients might appear to be at odds with the fact that population size histories can be inferred from allele frequency time series data ([@msw173-B40]; [@msw173-B2]; [@msw173-B16]). However, this is not the case. Methods for inferring the population size capitalize on information in the short-term fluctuations of the allele frequency around its expected value, arising from drift; conversely, estimators of selection coefficients capitalize on the long-term changes in allele frequency due to selection, effectively averaging over the short-term fluctuations due to drift. Our results suggest that allele frequencies often change quasi-deterministically, even in small populations. Thus, deviations around the expected trajectory can be distinguished from long-term changes, allowing effective population sizes to be inferred accurately even in small populations.

We have conducted our analyses under two different models of evolution: the discrete Wright--Fisher model and the continuous diffusion model. Although the diffusion model was developed as an approximation to the Wright--Fisher process, it also captures the limiting behavior of a large class of evolutionary models, including the Wright--Fisher process, as the population size grows to infinity and mutation and selection parameters are scaled accordingly. Thus, it is reasonable to believe that our findings will generalize to maximum likelihood estimators derived under a wide range of evolutionary models.

Taken together, our results help to characterize the properties of maximum likelihood methods for inferring selection coefficients from time series data. Because of the accuracy and beneficial properties of maximum likelihood methods, it is reasonable to believe that our results provide insight into the accuracy with which it is possible to infer selection coefficients from biological data, and the degree to which accounting for the true population history can improve these estimates. Our results also provide justification for the use of fast inference methods based on the assumption that allele frequencies evolve deterministically. Such methods can be applied to infer selection coefficients efficiently on large genomic data sets with many loci. Finally, our results provide further justification for the use of deterministic approximations in the development of statistical approaches for studying time series data.

Methods
=======

In this section, we compute the exact probability of an allele frequency trajectory in a population of piecewise-constant size under the discrete Wright--Fisher model and under the diffusion approximation. We also describe how drift can be ignored in these probabilities, yielding approximate estimators of selection coefficients that are similar to commonly used approaches that assume deterministic allele frequency trajectories.

Computing $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ under the Discrete Wright--Fisher Model
---------------------------------------------------------------------------------------------------------------------------------------------------

To compute the probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ under the discrete Wright--Fisher model, we make use of a hidden Markov model (HMM) similar to that presented in [@msw173-B48]. However, the hidden state in our discrete model is the count *c~t~* of the number of (unobserved) copies of allele *A* in the population at time *t*, rather than the continuous allele frequency *y~t~*.

In our model, the count *c~t~* of allele *A* evolves according to a Wright--Fisher process in which mutation occurs, followed by random mating, selection, and drift. Given that the count of allele *A* in generation *t* is *c~t ~*=*~ ~i*, let $f_{A|i}^{t}$ be the frequency of allele *A* in the gamete pool after mutation. Then $$f_{A|i}^{t} = \left( \frac{i}{2N_{t}} \right)(1 - u_{Aa}^{(t)}) + \left( {1 - \frac{i}{2N_{t}}} \right)u_{aA}^{(t)} = u_{aA}^{(t)} + (1 - u_{Aa}^{(t)} - u_{aA}^{(t)})\left( \frac{i}{2N_{t}} \right).$$

After random mating, the fraction of zygotes with each of the genotypes *AA*, *Aa*, and *aa* is ${(f_{A|i}^{t})}^{2},\, 2f_{A|i}^{t}(1 - f_{A|i}^{t})$, and ${(1 - f_{A|i}^{t})}^{2}$, from which it follows that the fraction of genotypes of each kind remaining after selection is given by $$\begin{array}{l}
{p_{AA|i}^{t} = \frac{{(f_{A|i}^{t})}^{2}(1 + s_{t})}{{\overline{w}}_{t,i}},} \\
{p_{Aa|i}^{t} = \frac{2f_{A|i}^{t}(1 - f_{A|i}^{t})(1 + h_{t}s_{t})}{{\overline{w}}_{t,i}},} \\
{p_{aa|i}^{t} = \frac{{(1 - f_{A|i}^{t})}^{2}}{{\overline{w}}_{t,i}},} \\
\end{array}$$ where ${\overline{w}}_{t,i} = {(f_{A|i}^{t})}^{2}(1 + s_{t}) + 2f_{A|i}^{t}(1 - f_{A|i}^{t})(1 + h_{t}s_{t}) + {(1 - f_{A|i}^{t})}^{2}$ is the mean fitness of the population.

Immediately after selection and before drift occurs, the probability that a randomly chosen allele is of type *A* is given by $p_{A|i}^{t} = p_{AA|i}^{t} + \frac{1}{2}p_{Aa|i}^{t}$. Then, as the result of drift, the count of allele *A* in generation *t* + 1 is binomially distributed with mean $p_{A|i}^{t}$. Thus, the probability that allele *A* has count *j* in generation *t* + 1, given that it had count *i* in generation *t* is $$\mathbb{P}_{\Theta,\mathcal{W}}\left\{ C_{t + 1} = j \middle| C_{t} = i \right\} = \left( \begin{array}{l}
{2N_{t + 1}} \\
\mathit{j} \\
\end{array} \right)\left( p_{A|i}^{t} \right)^{j}\left( 1 - p_{A|i}^{t} \right)^{2N_{t + 1} - j}.$$

The Wright--Fisher transition matrix $\mathbf{T}_{t,t + 1}$ from generation *t* to generation *t* + 1 is the $(2N_{t} + 1) \times (2N_{t + 1} + 1)$ matrix with entry *i*, *j* given by $${\lbrack\mathbf{T}_{t,t + 1}\rbrack}_{ij} = \mathbb{P}_{\Theta,\mathcal{W}}\left\{ C_{t + 1} = j \middle| C_{t} = i \right\},$$ which can be used to obtain the allele frequency distribution at each discrete generation *t* given the initial distribution at some time *r *\<* t*. In particular, define $\mathbf{d}_{t} = (\mathbb{P}\left\{ c_{t} = 0 \right\},\mathbb{P}\left\{ c_{t} = 1 \right\},\ldots,\mathbb{P}\left\{ c_{t} = 2N_{t} \right\})$, to be the distribution of the count of allele *A* in generation *t*. Using [equation (5)](#E13){ref-type="disp-formula"}, $\mathbf{d}_{t}$ can be computed recursively as $$\mathbf{d}_{t} = \mathbf{d}_{r}\left\lbrack {\prod\limits_{g = r + 1}^{t}\mathbf{T}_{g - 1,g}} \right\rbrack$$ for *r *\<* t*.

### Computing the Probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$

The probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ of the observed data is computed using the forward procedure for hidden Markov models. In particular, we define the vector $\mathbf{v}_{k}$ whose *i*th entry $v_{k,i}$ is the joint probability of the population-wide count of allele *A* at the *k*th sampling event and the observed sample allele counts up to sample *k*: $$v_{k,i} = \mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:k\rbrack} = o_{\lbrack 1:k\rbrack},C_{t_{k}} = i \right\}.$$ To simplify calculations, we also define the conditional "emission probability" $$\gamma_{i}(o_{k}) = \mathbb{P}_{\Theta}\left\{ O_{k} = o_{k} \middle| C_{t_{k}} = i \right\} = \left( \begin{array}{l}
n_{k} \\
o_{k} \\
\end{array} \right){(i/2N_{t_{k}})}^{o_{k}}{(1 - i/2N_{t_{k}})}^{n_{k} - o_{k}}$$ of the observed allele count, conditional on the population allele count. The probability in [equation (8)](#E16){ref-type="disp-formula"} comes from the fact that the observed allele count at time *t~k~* can be modeled as a binomial random variable with sample size *n~k~* and probability $c_{t_{k}}$. Although the observed allele count is, strictly speaking, hypergeometric we use the binomial distribution to maintain consistency with the formulas for the diffusion model, which are often derived using a binomial sampling distribution. The binomial and hypergeometric distributions are very similar for the population and sample sizes we consider. We then construct the emission probability vector $$\gamma(o_{k}) = (\gamma_{0}(o_{k}),\gamma_{1}(o_{k}),\ldots,\gamma_{2N_{t_{k}}}(o_{k})).$$

The probability of the data is then given by the forward procedure ([@msw173-B41]), outlined in Procedure 1. In Procedure 1, the formula for $\mathbf{v}_{1}$ comes from the fact that $$\begin{array}{l}
{\mathbf{v}_{1} = (\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{1} = o_{1},C_{t_{1}} = 0 \right\},\ldots,\,} \\
{\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{1} = o_{1},C_{t_{1}} = 2N_{\ell_{t_{1}}} \right\}) = (\gamma_{0}(o_{1})\mathbb{P}_{\Theta,\mathcal{W}}\left\{ C_{t_{1}} = 0 \right\},\ldots,} \\
{\gamma_{2N_{\ell_{t_{1}}}}(o_{1})\mathbb{P}_{\Theta,\mathcal{W}}\left\{ C_{t_{1}} = 2N_{\ell_{t_{1}}} \right\}) = \mathbf{d}_{t_{1}}\text{diag}\left\{ \gamma(o_{1}) \right\}\,} \\
{= \mathbf{d}_{0}\left\lbrack {\prod\limits_{t = 1}^{t_{1}}\mathbf{T}_{t - 1,t}} \right\rbrack\text{diag}\left\{ \gamma(o_{1}) \right\},} \\
\end{array}$$ where $\text{diag}(\gamma)$ denotes the square matrix whose diagonal entries are given by $\gamma$.

It has been noted by several authors that computing powers of the transition matrix is computationally prohibitive, providing one motivating factor for the use of approximations of the Wright--Fisher process, such as the diffusion and Gaussian approximations ([@msw173-B12]; [@msw173-B14]; [@msw173-B33]). However, the products $\prod_{t = t_{k - 1} + 1}^{t_{k}}\mathbf{T}_{t - 1,t}$ in Procedure 1 do not require products of the transition matrix $\mathbf{T}_{t - 1,t}$ because it suffices to repeatedly compute vector--matrix products instead of multiplying full matrices together. In practice, this can be done very quickly, even for large population sizes. A similar fast procedure was carried out by [@msw173-B56] to simulate trajectories under the Wright--Fisher model.

### Computing $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ When $\mathbf{d}_{0}$ Is Unspecified

When the initial distribution $\mathbf{d}_{0}$ is unspecified, the probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ can be obtained by assuming that the distribution $\mathbf{d}_{t_{1}}$ is uniform at the time of the first sampling event. Under this assumption, it follows directly from the second to last equality in [equation (10)](#E18){ref-type="disp-formula"} that the value of the joint density vector $\mathbf{v}_{1}$ is given by $$\mathbf{v}_{1} = \frac{1}{2N_{t_{1}}}1\text{diag}\left\{ \gamma(o_{1}) \right\} = \frac{1}{2N_{t_{1}}}\gamma(o_{1}),$$ where 1 is the vector of length $2N_{t_{1}}$ with all entries equal to one. This form of $\mathbf{v}_{1}$ can then be substituted into Procedure 1.

Computing $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ under the Diffusion Approximation
---------------------------------------------------------------------------------------------------------------------------------------------

The diffusion approximation models the evolution of the continuous population frequency *Y~t~* of allele *A*, rather than its count *C~t~*. The time-evolution of the random frequency *Y~t~* is governed by the diffusion transition density $p_{\Theta}(s,t;x,y)$ given by $$p_{\Theta}(s,t;x,y)dy = \,\mathbb{P}_{\Theta,\mathcal{D}}\left\{ y \leq Y_{t} < y + dy \middle| Y_{s} = x \right\},$$ for an infinitesimal increment *dy*. The quantity $p_{\Theta}(s,t;x,y)$ specifies the density of the allele frequency at time *t*, conditional on the value of the allele frequency at an earlier time *s*. For more details about the transition density function of the diffusion approximation, see Appendix A.

Using the diffusion transition density $p_{\Theta}(s,t;x,y)$ [@msw173-B48] developed an HMM to compute the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ of the data in a single epoch of constant size by efficiently integrating over the hidden allele frequencies $\left\{ y_{t_{1}},\ldots,y_{t_{K}} \right\}$ at the set of sampling times. Here, we extend this HMM to the case of piecewise-constant population size.

To compute the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ efficiently, [@msw173-B48] define the quantities $f_{k}(y)$ and $g_{k}(y)$ satisfying $$f_{k}(y)dy: = \,\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:k\rbrack} = o_{\lbrack 1:k\rbrack},y \leq Y_{t_{k}} < y + dy \right\},$$ and $$g_{k}(y)dy: = \,\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:k - 1\rbrack} = o_{\lbrack 1:k - 1\rbrack},y \leq Y_{t_{k}} < y + dy \right\}$$ for an infinitesimal increment *dy*. The quantity $f_{k}(y)$ is the joint density of the allele frequency at time *t~k~* and the observed counts up to sampling event *k*. The quantity $g_{k}(y)$ is the joint density of the allele frequency at time *t~k~* and the observed counts up to sampling event *k* -- 1.

It follows from the definition of $f_{k}(y)$ that the probability of the data is given by $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = {\int_{y = 0}^{1}f_{K}}(y)dy.$$

The quantity $f_{K}(y)$ can be obtained efficiently by recursion using the relationships $$f_{k}(y) = g_{k}(y){(\begin{array}{l}
n_{k} \\
o_{k} \\
\end{array})}y^{o_{k}}{(1 - y)}^{n_{k} - o_{k}},$$ and $$g_{k}(y) = {\int_{z = 0}^{1}f_{k - 1}}(z)p_{\Theta}(t_{k - 1},t_{k};z,y)dz.$$

[Equation (16)](#E24){ref-type="disp-formula"} follows from the fact that the observed number of copies of allele *A* at sampling event *k* is binomially distributed with count *n~k~* and probability $y_{t_{k}}$ and [equation (17)](#E25){ref-type="disp-formula"} follows from the law of total probability integrating over $Y_{t_{k - 1}}$.

Let $B_{\ell,i}(y)$ $(i = 0,1,\ldots)$ be the *i*th eigenfunction of the backward diffusion operator $\ell_{\ell}$ and let $\pi_{\ell}(y)$ be the speed density of $\ell_{\ell}$ in Epoch $\ell$ (Appendix A). [@msw173-B48] demonstrated that the recursive formulas in [equations (16)](#E24){ref-type="disp-formula"} and [(17)](#E25){ref-type="disp-formula"} can be evaluated efficiently by expressing $f_{k}(y)$ and $g_{k}(y)$ as series of the form $$f_{k}(y) = \sum\limits_{i = 0}^{\infty}b_{k,i}\pi_{\ell_{k}}(y)B_{\ell_{k},i}(y) = \mathbf{b}_{k}\pi_{\ell_{k}}(y)\mathbf{B}_{\ell_{k}}(y)$$ and $$g_{k}(y) = \sum\limits_{i = 0}^{\infty}a_{k,i}\pi_{\ell_{k}}(y)B_{\ell_{k},i}(y) = \mathbf{a}_{k}\pi_{\ell_{k}}(y)\mathbf{B}_{\ell_{k}}(y),$$ where $\mathbf{B}_{\ell}(y) = (B_{\ell,0}(y),\, B_{\ell,1}(y),\,\ldots)$, and where $\mathbf{b}_{k} = (b_{k,0},b_{k,1},\ldots)$ and $\mathbf{a}_{k} = (a_{k,0},a_{k,1},\ldots)$ are vectors of constants that encode the densities $f_{k}(y)$ and $g_{k}(y)$. In Appendix B, we extend the results of [@msw173-B48] to derive recursive formulas for the coefficients $\mathbf{a}_{k}$ and $\mathbf{b}_{k}$, resulting in Procedure 2, which computes the probability of an allele frequency trajectory under the diffusion approximation in a population of piecewise constant size.

### Computing $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ When y~0~ Is Unspecified

When the initial frequency *y*~0~ is unspecified, the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ can be obtained by assuming that the distribution of the allele frequency $y_{t_{1}}$ at the time of the first sampling event is uniform. Under this assumption, we show in Lemma B.3.2 that the coefficients $\mathbf{b}_{1}$ encoding the distribution at the time of the first sampling event are given by [equation (B.17)](#E82){ref-type="disp-formula"}. The form of $\mathbf{b}_{1}$ in [equation (B.17)](#E82){ref-type="disp-formula"} is then used in Procedure 2.

Conditional Probabilities
-------------------------

Sometimes it is desirable to compute the probability of the observed allele counts conditional on the event *S~K~* that allele *A* is segregating in the final sample. In this section, we provide formulas for these conditional probabilities under the Wright--Fisher and diffusion models.

### Computing $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}$

In the "Simulations Conditioning on Segregation" section, we consider the probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}$ of the data conditional on the event *S~K~* that allele *A* is segregating in the final sample. In Appendix C, we show that in the case of the discrete Wright--Fisher model, $$\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}\, = \frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \right\}}\sum\limits_{i = 0}^{2N_{t_{K}}}v_{K,i},$$ where $v_{k,i}$ is defined in [equation (7)](#E15){ref-type="disp-formula"} and $\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\} = 1$ if $1 \leq o_{K} < n_{K}$, or 0 otherwise. The probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \right\}$ is given in [equation (C.3)](#E94){ref-type="disp-formula"}. Thus, if we wish to compute conditional probabilities under the Wright--Fisher model, we carry out Procedure 1, replacing step 3 with [equation (20)](#E28){ref-type="disp-formula"}.

### Computing $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| F_{K} \right\}$

Similarly, for the event *F~K~* that allele *A* is segregating or fixed in the final sample, we show in Appendix C that $$\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| F_{K} \right\}\, = \frac{\mathbb{P}\left\{ F_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{W}}\left\{ F_{K} \right\}}\sum\limits_{i = 0}^{2N_{t_{K}}}v_{K,i},$$ where $v_{k,i}$ is defined in [equation (7)](#E15){ref-type="disp-formula"} and $\mathbb{P}\left\{ F_{K} \middle| O_{K} = o_{K} \right\} = 1$ if $1 \leq o_{K} \leq n_{K}$, or 0 otherwise. The probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ F_{K} \right\}$ is given in [equation (C.6)](#E97){ref-type="disp-formula"}. If we wish to compute the conditional probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| F_{K} \right\}$ under the Wright--Fisher model, we carry out Procedure 1, replacing step 3 with [equation (21)](#E29){ref-type="disp-formula"}.

### Computing $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}$

In the case of the diffusion approximation, we show in Appendix D that the conditional probability of the data given *S~K~* can be computed as $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}\, = \frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\} c_{\ell_{K},0}b_{K,0}}{B_{\ell_{K},0}\left( 0 \right) - c_{\ell_{K},0}{\widetilde{b}}_{K,0}\left( 0 \right) - c_{\ell_{K},0}{\widetilde{b}}_{K,0}\left( n_{K} \right)},$$ where $\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\} = 1$ if $1 \leq o_{K} < n_{K}$ or 0 otherwise, and $${\widetilde{\mathbf{b}}}_{K}\left( j \right) = \left\{ \begin{array}{l}
{\mathbf{b}_{0}\mathbf{E}_{\ell_{1}}\left( t_{K} \right)\mathbf{W}_{\ell_{K}}\mathbf{G}_{\ell_{K}}^{j}\left( 1 - \mathbf{G}_{\ell_{K}} \right)^{n_{K} - j}\mathbf{W}_{\ell_{K}}^{- 1},} \\
{\,\text{if}\,\ell_{t_{K}} = 1,} \\
{\mathbf{b}_{0}\mathbf{F}\left( 0,t_{K};\zeta \right)\mathbf{W}_{\ell_{K}}\mathbf{G}_{\ell_{K}}^{j}\left( 1 - \mathbf{G}_{\ell_{K}} \right)^{n_{K} - j}\mathbf{W}_{\ell_{K}}^{- 1},} \\
{\,\text{otherwise}\,.} \\
\end{array} \right.\ $$

Thus, if we are interested in conditional probabilities under the diffusion model, we carry out Procedure 2, replacing step 3 with [equation (22)](#E30){ref-type="disp-formula"}.

The Probability in the Absence of Genetic Drift
-----------------------------------------------

If we ignore genetic drift, the allele frequency changes deterministically over time, as it would in a population of infinite size. Here, we obtain versions of Procedures 1 and 2 in the case when the changes in allele frequency arising from genetic drift are negligible relative to the changes due to selection and recurrent mutation.

### Deterministic Allele Frequency Trajectories under the Wright--Fisher Model

If there is no contribution to the change in allele frequency arising from genetic drift, the allele frequency in a given generation is equal to its expectation after mutation, random mating, and selection, conditional on its value in the previous generation. Because the expectation is not necessarily integer-valued, we no longer consider discrete integer allele counts *c~t~*. Instead, we track the expected allele frequency in the absence of drift, which we denote by $y_{t}^{\infty} \equiv \mathbb{E}_{\infty}\lbrack Y_{t}\rbrack$, where the subscript ∞ denotes the expectation without drift as the effective population size tends to infinity.

The expected frequency $y_{t}^{\infty}$ in the absence of drift is obtained by combining [equations (2)](#E10){ref-type="disp-formula"} and [(3)](#E11){ref-type="disp-formula"}, ignoring the drift step in [equation (4)](#E12){ref-type="disp-formula"}, yielding $$y_{t + 1}^{\infty} = \left\lbrack \frac{{({\widetilde{y}}_{t}^{\infty})}^{2}(1 + s_{t}) + {\widetilde{y}}_{t}^{\infty}(1 - {\widetilde{y}}_{t}^{\infty})(1 + h_{t}s_{t})}{{\overline{w}}_{t}} \right\rbrack,$$ where $${\widetilde{y}}_{t}^{\infty} = u_{aA}^{(t)} + (1 - u_{Aa}^{(t)} - u_{aA}^{(t)})y_{t}^{\infty}$$ and ${\overline{w}}_{t} = {({\widetilde{y}}_{t}^{\infty})}^{2}(1 + s_{t}) + 2{\widetilde{y}}_{t}^{\infty}(1 - {\widetilde{y}}_{t}^{\infty})(1 + h_{t}s_{t}) + {(1 - {\widetilde{y}}_{t}^{\infty})}^{2}$. [Equations (24)](#E32){ref-type="disp-formula"} and [(25)](#E33){ref-type="disp-formula"} are iterated to find the allele frequency in any generation *t* \> 0.

[Equation (24)](#E32){ref-type="disp-formula"} is related to [equation 3.1](#E11){ref-type="disp-formula"} of [@msw173-B21], which describes the dynamics of allele frequency change in a population of infinite size without drift. The frequency trajectories described by this formula are closely approximated by logistic curves that have closed formulas ([@msw173-B14], [equation 3](#E11){ref-type="disp-formula"}), which can be used to increase the speed of deterministic approaches even further. However, we have chosen to implement the exact formulas for better comparison with the exact likelihoods.

### Deterministic Allele Frequency Trajectories under the Diffusion Model

Under the diffusion model in an Epoch $\ell$ of constant size, the allele frequency *Y~t~* obeys the stochastic differential equation (SDE) $$dY_{t} = \mathcal{M}_{\ell}(Y_{t})dt + \sqrt{Y_{t}(1 - Y_{t})}dB_{t},\, t \in \lbrack\tau_{\ell - 1},\tau_{\ell}\rbrack,$$ with the initial condition $Y_{\tau_{\ell - 1}} = y_{\tau_{\ell - 1}}$, where time is measured in units of generations and $\tau_{\ell - 1}$ is the time at which Epoch $\ell$ begins ([@msw173-B31], Section 15.14; [@msw173-B11], Section 7.2). The quantity $\sqrt{Y_{t}(1 - Y_{t})}$ in [equation (26)](#E34){ref-type="disp-formula"} controls random fluctuations due to drift whereas the quantity $\mathcal{M}_{\ell}(y)$ describes the deterministic change in the mean frequency of the allele over time due to mutation and selection and is given by $$\mathcal{M}_{\ell}(y) = u_{aA}^{(\ell)} - (u_{aA}^{(\ell)} + u_{Aa}^{(\ell)})y + \, y(1 - y)\lbrack h_{\ell}s_{\ell}(1 - 2y) + s_{\ell}y\rbrack.$$ In [equation (27)](#E35){ref-type="disp-formula"}, we have rescaled the usual form of $\mathcal{M}_{\ell}$ so that time is measured continuously in units of generations.

If the drift term in [equation (26)](#E34){ref-type="disp-formula"} is negligible compared with $\mathcal{M}_{\ell}(Y_{t})$, then [equation (26)](#E34){ref-type="disp-formula"} can be approximated by the ordinary differential equation $$\frac{dy_{t}^{\infty}}{dt} = \mathcal{M}_{\ell}(y_{t}^{\infty}),$$ where we may write $y_{t}^{\infty}$ instead of *Y~t~* because the evolution of the allele frequency is deterministic and follows its expectation in the absence of drift.

We can also suppress the explicit dependence on the epoch $\ell$ by defining $\mathcal{M}_{t}(y_{t}^{\infty}) \equiv \mathcal{M}_{\ell_{t}}(y_{t}^{\infty})$, yielding $$\frac{dy_{t}^{\infty}}{dt} = \mathcal{M}_{t}(y_{t}^{\infty}),y_{0}^{\infty} = y_{0},\,\, t \in \lbrack 0,\tau_{L}\rbrack,$$ which holds for the full population history across all epochs $\ell = 1,\ldots,L$. [Equation (29)](#E37){ref-type="disp-formula"} can be solved numerically, for instance by choosing a sufficiently small time step $\bigtriangleup t$ and iteratively computing $y_{t + \bigtriangleup t}^{\infty} = \mathcal{M}_{t}(y_{t}^{\infty}) \bigtriangleup t$.

### Sample Probabilities Based on Deterministic Allele Frequency Trajectories

To compute the probability $\mathbb{P}_{\Theta}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ under either the discrete Wright--Fisher or diffusion models when drift is negligible, we note that the observations $(O_{1},\ldots,O_{K})$ are conditionally independent of one another, given the underlying allele frequencies. Thus, in the absence of drift we have $$\mathbb{P}_{\Theta}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}\, = \prod\limits_{k = 1}^{K}\mathbb{P}_{\Theta}\left\{ O_{k} = o_{k} \middle| Y_{t_{k}} = y_{t_{k}}^{\infty} \right\}$$ for both the diffusion and Wright--Fisher models, where $y_{t_{k}}^{\infty}$ is the deterministic allele frequency at time *t~k~*, for $k = 1,\ldots,K$. Using [equations (24)](#E32){ref-type="disp-formula"} and [(30)](#E38){ref-type="disp-formula"}, the probability of the data under the Wright--Fisher model in a population without drift can be obtained using Procedure 3. Similarly, using [equations (29)](#E37){ref-type="disp-formula"} and [(30)](#E38){ref-type="disp-formula"}, the probability of the data in the case of the diffusion model is given by Procedure 4.

### Deterministic Estimators When the Initial Allele Frequency Is Unknown

If the initial allele frequency is unspecified, the selection coefficient can still be inferred using the deterministic estimators by integrating over the allele frequency $y_{t_{1}}$ at the time of the first sampling event, where we assume in these analyses that the initial allele frequency distribution is uniform. In practice, we perform this integration by specifying a dense grid of values $y_{t_{1}} \in \lbrack 0, \bigtriangleup y,2 \bigtriangleup y,\ldots,1\rbrack$, where $\bigtriangleup y = 1/M$ for some large predetermined value *M*. Defining the probability of the data when the initial allele frequency is *y* to be $\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack};y_{0} = y \right\}$, we compute the probability of the data as $$\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}\, = \frac{1}{M + 1}\sum\limits_{j = 0}^{M}\mathbb{P}_{\Theta,\mathcal{W}}^{\infty}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack};y_{0} = j \bigtriangleup y \right\}.$$

Simulations
-----------

Allele frequency trajectories were simulated under two different models, the discrete Wright--Fisher model and the continuous diffusion model. All simulations were carried out by iteratively sampling the allele frequency at the time points $t_{1},t_{2},\ldots,t_{K}$ starting with a specified frequency *y*~0~ or with the initial frequency *y*~0~ sampled uniformly from the interval $(0,0.5\rbrack$ at time $t_{0} = 0$.

Wright--Fisher simulations were carried out using Procedure 5. Specifically, for an initial allele frequency *y*~0~, we took the initial distribution $\mathbf{d}_{0}$ to be the standard basis vector $\mathbf{e}_{\lceil y_{0}\rceil}$ of length $2N_{0} + 1$ with element $\lceil{2N_{0}y_{0}}\rceil$ set to unity and all other elements set to zero. We then sampled the population allele count $c_{t_{k}}$ at each sampling time *t~k~* by iteratively propagating the allele frequency distribution $d_{t_{k}}$ forward, conditional on the population count $c_{t_{k - 1}}$ at the previous sampling time using [equation (6)](#E14){ref-type="disp-formula"}. The derived allele count *o~k~* in each sample *k* was then chosen from a binomial distribution with sample size *n~k~* and probability $c_{t_{k}}/2N_{t_{k}}$. For the results in the "Inference Accuracy for Populations of Constant Finite Size", "Inference Accuracy in Populations of Piecewise Constant Size", "Conditioning on Segregation in the Final Sample", and "The Effect of Sample Size on Accuracy" sections, we fixed $y_{0} = 0.01$ or $y_{0} = 0.1$ in all simulations, as indicated in the results section. In the "Unspecified Initial Allele Frequencies" section, we simulated *y*~0~ uniformly from the interval $(0,0.5\rbrack$.

Simulations under the diffusion model were carried out using Procedure 6. Specifically, for an initial frequency *y*~0~, we computed the coefficients $\mathbf{b}_{0}$ of the expansion of the initial condition $\delta(y - y_{0})$ in the basis functions $\mathbf{B}_{\ell_{1}}$ of the first epoch. In Appendix D, we show that the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{k} = i \right\}$ of sampling *i* derived alleles at sampling event *k* (eq. D.3) can be obtained by computing a set of coefficients that we call ${\widetilde{\mathbf{b}}}_{k}(i)$ (eq. D.7) that specify an expansion in the basis of eigenfunctions $\mathbf{B}_{\ell_{k}}$.

Thus, to sample the derived allele count *o*~1~ in the first sample, we propagated the initial coefficients $\mathbf{b}_{0}$ to obtain ${\widetilde{\mathbf{b}}}_{1}(i)$ and we used the fact that $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{1} = i \right\} \propto {\widetilde{\mathbf{b}}}_{1,0}(i)$ (eq. D.3) to sample $o_{1} = i$ with probability ${\widetilde{\mathbf{b}}}_{1}(i)/\sum\limits_{i = 0}^{n_{1}}{\widetilde{\mathbf{b}}}_{1}(i)$. Using the fact that ${\widetilde{\mathbf{b}}}_{k}(o_{k})$ is defined in Section D as the set of coefficients such that $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{k} = o_{k},y \leq Y_{t_{k}} < y + dy \right\} = {\widetilde{\mathbf{b}}}_{k}(o_{k})\pi_{\ell_{k}}\mathbf{B}_{\ell_{k}}(y)$ for an infinitesimal increment *dy*, we then set $\mathbf{b}_{1} = {\widetilde{\mathbf{b}}}_{1}(o_{k})$ and iterated this procedure to obtain the subsequent coefficient vectors $\left\{ {\widetilde{\mathbf{b}}}_{k}(i) \right\}_{i = 2}^{n_{k}}$ and the full set of samples $o_{\lbrack 1:K\rbrack}$. Simulated trajectories were checked for accuracy against those obtained using the software of [@msw173-B29], personal communication.

Grid Searches to Infer Selection Coefficients
---------------------------------------------

The same iterated grid search procedure, implemented as a wrapper function and outlined in Procedure 7, was used to infer selection coefficients using each of the methods described in Procedures 1--4. The grid search is initialized by specifying lower and upper bounds on a region $\lbrack s_{\ell},s_{u}\rbrack$ over which to search for a value of the selection coefficient *s* that maximizes the likelihood. Assuming that the likelihood surface is convex and smooth, the grid search iteratively refines the search to the region that contains the optimum until the width of the region is smaller than a specified tolerance *ϵ*. In our analyses, we chose the initial region to be $\lbrack s_{\ell},s_{u}\rbrack = \lbrack - 0.99,1\rbrack$ and we set $\epsilon = 2 \cdot 10^{- 4}$.

For each allele frequency trajectory simulated under the Wright--Fisher model, Procedure 7 was carried out using Procedure 1 and then subsequently for the same trajectory using Procedure 3, with or without conditioning. For each allele frequency trajectory simulated under the diffusion model, Procedure 7 was carried out using Procedure 2 and then subsequently for the same trajectory using Procedure 4, with or without conditioning. A gradient descent optimization approach produced nearly identical results and is available as an option in the software package we have released.

Computing the Watterson Estimator of *N~e~* from the Expected SFS of a Piecewise Constant Population
----------------------------------------------------------------------------------------------------

In the "Inference Accuracy in Populations of Piecewise Constant Size" section, we obtain a crude constant estimate of the effective size *N~e~* of a piecewise constant population by computing the expected unnormalized SFS for the true history and then inferring the effective size of a constant population with the same level of genetic diversity as the piecewise constant population using the Watterson estimator.

Under the assumption that any given base pair in a collection of *n* haplotypes from a population is at most biallelic, the unnormalized and folded SFS for the *n* haplotypes is a collection of counts $\left\{ \xi_{n,i} \right\}_{i = 1}^{n - 1}$ in which $\xi_{n,i}$ is the number of sites with one or the other allele appearing in *i* out of *n* copies in the sample. When the expected number *μ* of new mutations occurring in a population within in a region of fixed length is specified, the expected SFS for *n* haplotypes spanning the region can be computed for a population of piecewise constant size using algorithms by [@msw173-B30].

The effective size of a population of constant size extending infinitely far back into the past with the same level of diversity as the piecewise constant population can be computed using the Watterson estimator ([@msw173-B54], Eqn. 1.4a; [@msw173-B27], p. 62). The information used by the Watterson estimator is the diversity in the *n* sampled haplotypes of specified length, combined with knowledge of the expected mutation rate *μ*. The Watterson estimator is given by $${\widehat{N}}_{e} = \frac{1}{4\mu}\frac{S}{\sum\limits_{i = 1}^{n - 1}i},$$ where *S* is the total number of segregating sites observed in the sample. The total number of segregating sites can be computed as $S = \sum_{i = 1}^{n - 1}\xi_{n,i}$, where $\xi_{n,i}$ is the *i*th entry of the un-normalized folded SFS.

Thus, given the expected SFS $\left\{ \xi_{n,i} \right\}_{i = 1}^{n - 1}$ for a piecewise constant population computed using the algorithm of [@msw173-B30] and a specified mutation rate *μ*, we computed the Watterson estimate of the size *N~e~* of a constant-sized population with similar diversity as $${\widehat{N}}_{e} = \frac{1}{4\mu}\frac{\sum\limits_{i = 1}^{n - 1}\xi_{n,i}}{\sum\limits_{i = 1}^{n - 1}i}.$$

When computing the expected SFS and the Watterson estimator, we arbitrarily chose a value of *μ*  = 1. **Procedure 5.** Sampling Wright--Fisher trajectories1: Optional: sample the initial population allele frequency y~0~ uniformly from a predetermined interval $\left( y_{\ell},y_{u} \right\rbrack$, where $0 \leq y_{\ell} \leq y_{u} \leq 1$, and set $c_{t_{0}} = \left\lceil {2N_{0}y_{0}} \right\rceil$.2: Let $\mathbf{d}_{0} = \mathbf{e}_{c_{t_{0}}}$ be the standard basis vector of length $2N_{0} + 1$ with element $c_{t_{0}}$ equal to one.3: For $k = 1:K$,Compute the conditional distribution ${\widetilde{\mathbf{d}}}_{k}$ at sampling event *k*, conditional on the distribution $\mathbf{d}_{k - 1}$ at time $t_{k - 1}$ using $${\widetilde{\mathbf{d}}}_{k} = \mathbf{d}_{k - 1}\left\lbrack {\prod\limits_{t = t_{k - 1} + 1}^{t_{k}}\mathbf{T}_{t - 1,t}} \right\rbrack.$$Sample $c_{t_{k}} = i$ with probability ${\widetilde{\mathbf{d}}}_{k}\left( i \right)$, for $i = 0,1,\ldots,2N_{t_{k}}$.Set $\mathbf{d}_{k} = \mathbf{e}_{c_{t_{k}}}$, where $\mathbf{e}_{c_{t_{k}}}$ is the standard basis vector of length $2N_{t_{k}} + 1$ with element $c_{t_{k}}$ equal to one.4: For$k = 1:K$, sample $O_{k} \sim \text{Binomial}\left( n_{k},c_{t_{k}}/\left( 2N_{t_{k}} \right) \right)$.5: If conditioning on segregation in the final sample and *o~K~* = 0 or *o~K ~*=*~ ~n~K~*, return to Step 1. Otherwise, the sampled trajectory is $o_{\lbrack 1:K\rbrack}$. **Procedure 6.** Simulating diffusion model trajectories1: Optional: sample the initial population allele frequency y~0~ uniformly from a predetermined interval $\left( y_{\ell},y_{u} \right\rbrack$, where $0 \leq y_{\ell} \leq y_{u} \leq 1$.2: For an initial starting frequency y~0~ initialize $$\mathbf{b}_{0} = \mathbf{C}_{\ell_{1}}^{- 1}\mathbf{B}_{\ell_{1}}\left( y_{0} \right),$$ where $\mathbf{B}_{\ell}\left( y_{0} \right)$ is the vector of eigenfunctions of the diffusion operator given inequation (A.14) and $\mathbf{C}_{\ell} = \text{diag}\left\{ \left\langle B_{\ell,i},B_{\ell,i} \right\rangle \right\}_{i = 0}^{\infty}$ is given in equation(A.18).3: For $k = 1:K$,Compute $$\mathbf{a}_{k}~ = \begin{cases}
{\mathbf{b}_{k - 1}\mathbf{E}_{\ell_{k}}\left( t_{k} - t_{k - 1} \right)} & {\,\text{if}\,\ell_{k - 1} = \ell_{k},} \\
{\mathbf{b}_{k - 1}\mathbf{F}\left( t_{k - 1},t_{k};\zeta \right)} & {\,\text{otherwise}.} \\
\end{cases}~~$$For $i = 0,1,\ldots,n_{k}$, compute $${\widetilde{\mathbf{b}}}_{k}\left( i \right) = \mathbf{a}_{k}\mathbf{W}_{\ell_{k}}\mathbf{G}_{\ell_{k}}^{i}\left( 1 - \mathbf{G}_{\ell_{k}} \right)^{n_{k} - i}\mathbf{W}_{\ell_{k}}^{- 1},$$ noting that ${\widetilde{\mathbf{b}}}_{k,0}\left( i \right) \propto \mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{k} = i \right\}$ (eq. D.3).Sample o~*k*~ = i with probability $\mathbf{b}_{k,0}\left( i \right)/\sum\limits_{i = 0}^{n_{k}}\mathbf{b}_{k,0}\left( i \right)$.Set $\mathbf{b}_{k} = \mathbf{b}_{k}\left( o_{k} \right)$.Here, the matrices $\mathbf{E}_{\ell}\left( t \right),\,\mathbf{F}\left( t_{k - 1},t_{k};\zeta \right),\,\mathbf{W}_{\ell}$, and $\mathbf{G}_{\ell}$ are given by equations (A.17), (B.10), (A.15) and (A.11), respectively and ζ is the set of Chebyshev nodes in the interval $\left\lbrack 0,1 \right\rbrack$. The matrix inverse $\mathbf{W}_{\ell}^{- 1} = \mathbf{D}_{\ell}\mathbf{W}_{\ell}^{T}\mathbf{C}_{\ell}^{- 1}$ is computed easily using the diagonal matrices $\mathbf{C}_{\ell}$ and $\mathbf{D}_{\ell}$ in equations (A.18) and (A.19).4: If conditioning on segregation in the final sample and *o~K~* = 0 or *o~K ~*=*~ ~n~K~*, return to Step 1. Otherwise, the sampled trajectory is $o_{\lbrack 1:K\rbrack}$. **Procedure 7.** Grid search1: Specify the bounds $s_{\ell}$ and s~u~ of an interval $\lbrack s_{\ell},s_{u}\rbrack$in which to search.2: Specify the stopping tolerance ε.3: While $s_{u} - s_{\ell} > \varepsilon$:Evaluate the likelihood at the points $\left\{ s^{(0)},s^{(1)},s^{(2)},s^{(3)},s^{(4)},s^{(5)} \right\} = \left\{ s_{\ell} + i\delta \right\}_{i = 0}^{5}$, where $\delta = (s_{u} - s_{\ell})/5$.Let $s_{\max\,} = \arg\,{\max\,}_{i}\left\{ s^{(i)} \right\}$.If $s_{\max\,} = 0$, set $s_{\ell} = s^{(0)}$ and $s_{u} = s^{(1)}$.Else, if $s_{\max\,} = 5$, set $s_{\ell} = s^{(4)}$ and $s_{u} = s^{(5)}$.Else, set $s_{\ell} = s^{(s_{\max\,} - 1)}$ and $s_{u} = s^{(s_{\max\,} + 1)}$.4: Return $\widehat{s} = (s_{\ell} + s_{u})/2$.

Supplementary Material
======================

[Supplementary figure S1](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msw173/-/DC1) is available at *Molecular Biology and Evolution* online (<http://www.mbe.oxfordjournals.org/>).

###### Supplementary Data

This research was supported by the National Institutes of Health (grant number R01-GM094402) and by a Packard Fellowship for Science and Engineering. We would also like to thank Jeffrey Spence, Joshua Schraiber, and an anonymous reviewer for their helpful comments on this manuscript.

Diffusion Transition Densities: Background
==========================================

The equations in the "Computing $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ under the Diffusion Approximation" section were derived under a model in which the selected allele *A* evolves under the diffusion approximation in a population of piecewise constant size. Given that allele *A* has frequency *x* at a fixed time *s*, the density at a later time *t* is given by the transition density of the diffusion approximation ([eq. 12](#E20){ref-type="disp-formula"}). [@msw173-B48] derived a formula for the density for the case of a single population of constant size. Here, we review this derivation to provide background and notation for the derivation of the diffusion model probability computed in Procedure 2.

The Diffusion Approximation in a Population of Constant Size
------------------------------------------------------------

Let $p_{\ell}(s,t;x,y)$ denote the transition density restricted to a specific epoch $\ell$ of constant size with $s,t \in \ell$. The density $p_{\ell}(s,t;x,y)$ is the unique solution of the Kolmogorov backward equation, $$\frac{\partial p_{\ell}(s,t;x,y)}{\partial t} = \frac{1}{2N_{\ell}}\mathcal{L}_{\ell}p_{\ell}(s,t;x,y)$$ satisfying the terminal condition $\rho_{s}(y) = \delta(y - x)$, where $\delta( \cdot )$ is the Dirac delta distribution and $\mathcal{L}_{\ell}$ is the Kolmogorov backward operator in the epoch defined in [equation (A.2)](#E47){ref-type="disp-formula"}. The factor $1/2N_{\ell}$ in [equation (A.1)](#E46){ref-type="disp-formula"} is introduced so that the time-scaling is the same in all epochs, and time is measured continuously in units of generations.

The Kolmogorov backward operator is defined in terms of the scaled mutation and selection parameters $\beta_{\ell} = 4N_{\ell}u_{aA}^{(\ell)},\,\alpha_{\ell} = 4N_{\ell}u_{Aa}^{(\ell)}$, and $\sigma_{\ell} = N_{\ell}s_{\ell}$ as $$\left. \mathcal{L}_{\ell} = \frac{1}{2}\xi^{2}(x)\frac{\partial^{2}}{\partial x^{2}} + \mu_{\ell}(x \right)\frac{\partial}{\partial x},$$ where the quantity $$\xi^{2}(x) = x(1 - x)$$ captures the contribution to the change in allele frequency arising from genetic drift and $$\mu_{\ell}(x) = \frac{1}{2}\left\lbrack \beta_{\ell} - (\beta_{\ell} + \alpha_{\ell})x\rbrack + \, 2x(1 - x)\lbrack h_{\ell}\sigma_{\ell}(1 - 2x) + \sigma_{\ell}x\rbrack \right.$$ captures the contribution from recurrent mutation and selection.

[@msw173-B47] showed that $p_{\ell}(s,t;x,y)$ can be expressed as an expansion in the eigenfunctions of $\mathcal{L}_{\ell}$ of the form $$p_{\ell}(s,t;x,y) = \,\sum\limits_{n = 0}^{\infty}e^{- \lambda_{\ell,n}(t - s)/2N_{\ell}}\frac{\pi_{\ell}(y)B_{\ell,n}(x)B_{\ell,n}(y)}{\left\langle B_{\ell,n},B_{\ell,n} \right\rangle_{\pi_{\ell}}},$$ where $\left\{ B_{\ell,n}(x) \right\}_{n = 0}^{\infty}$ are the eigenfunctions of $\ell_{\ell}$ with associated eigenvalues $\left\{ \lambda_{\ell,n} \right\}_{n = 0}^{\infty}$ and the function $\pi_{\ell}(y)$ is given by $$\pi_{\ell}(y) = e^{{\overline{\sigma}}_{\ell}(y)}y^{\beta_{\ell} - 1}{(1 - y)}^{\alpha_{\ell} - 1},$$ where ${\overline{\sigma}}_{\ell}(y) = 4h_{\ell}\sigma_{\ell}y(1 - y) + 2\sigma_{\ell}y^{2}$. The inner product $\left\langle f,g \right\rangle_{\omega}$ with respect to a weight function $\omega(x)$ in [equation (A.5)](#E50){ref-type="disp-formula"} is defined for two functions *f* and *g* on an interval $\lbrack a,b\rbrack$ by $$\left\langle f,g \right\rangle_{\omega} = {\int_{a}^{b}f}(x)g(x)\omega(x)dx.$$

In equation (A.5), the inner product $\left\langle \cdot , \cdot \right\rangle_{\pi_{\ell}}$ is taken over the interval $\lbrack 0,1\rbrack$ with respect to $\pi_{\ell}(y)$.

Expressions for the Quantities in Equation (A.5)
================================================

Expressions for the eigenvalues $\left\{ \lambda_{\ell,n} \right\}_{n = 0}^{\infty}$, eigenfunctions $\left\{ B_{\ell,n}(y) \right\}_{n = 0}^{\infty}$, and inner products $\left\{ \left\langle B_{\ell,n},B_{\ell,n} \right\rangle \right\}_{n = 0}^{\infty}$ in [equation (A.5)](#E50){ref-type="disp-formula"} can be obtained using a matrix formulation developed by [@msw173-B48]. In particular, the eigenfunctions $\left\{ B_{\ell,n}(y) \right\}_{n = 0}^{\infty}$ can be expressed as $$B_{\ell,n}(y) = \sum\limits_{m = 0}^{\infty}w_{\ell,n,m}e^{- {\overline{\sigma}}_{\ell}(y)/2}R_{m}^{(\beta_{\ell},\alpha_{\ell})}(y),$$ where $R_{m}^{(\alpha,\beta)}(y) = p_{m}^{(\beta - 1,\alpha - 1)}(2y - 1)$ and $p_{m}^{(a,b)}(y)$ is the *m*th classical Jacobi polynomial ([@msw173-B1], Chapter 22). The vector $\mathbf{w}_{\ell,n} = (w_{\ell,n,0},w_{\ell,n,1},\ldots)$ is the *n*th left eigenvector of the infinite-dimensional matrix $$\mathbf{M}_{\ell}: = - {({\Upsilon^{(\alpha_{\ell},\beta_{\ell})} + \sum\limits_{r = 0}^{4}q_{\ell,r}\mathbf{G}_{\ell}^{r}})}$$ corresponding to the *n*th eigenvalue $\lambda_{\ell,n}$, where $\Upsilon^{(\alpha,\beta)} = \text{diag}(\upsilon_{0}^{(\alpha,\beta)},\upsilon_{1}^{(\alpha,\beta)},\ldots)$ is the diagonal matrix with elements given by $\upsilon_{n}^{(\alpha,\beta)} = \frac{1}{2}n(n + \alpha + \beta - 1)$ and the quantities $q_{\ell,r}$ are given by $$\begin{array}{l}
{q_{\ell,0} = \alpha_{k}h_{\ell}\sigma_{\ell},} \\
{q_{\ell,1} = - (2 + 3\alpha_{\ell} + \beta_{\ell} - 2h_{\ell}\sigma_{\ell})h_{\ell}\sigma_{\ell} + (1 + \alpha_{\ell})\sigma_{\ell},} \\
{q_{\ell,2} = (2 + 2\alpha_{\ell} + 2\beta_{\ell} + 4\sigma_{\ell} - 10h_{\ell}\sigma_{\ell})h_{\ell}\sigma_{\ell} - \,(1 + \alpha_{\ell} + \beta_{\ell})\sigma_{\ell},} \\
{q_{\ell,3} = 16h_{\ell}^{2}\sigma_{\ell}^{2} + 2\sigma_{\ell}^{2}(1 - 6h_{\ell}),} \\
{q_{\ell,4} = - 2\sigma_{\ell}^{2}{(1 - 2h_{\ell})}^{2}.} \\
\end{array}$$

The matrix $\mathbf{G}_{\ell}^{r}$ in [equation (A.9)](#E54){ref-type="disp-formula"} has elements given by $${\lbrack\mathbf{G}_{\ell}\rbrack}_{n,m}\  = \left\{ \begin{array}{l}
{\frac{(n + \alpha_{\ell} - 1)(n + \beta_{\ell} - 1)}{(2n + \alpha_{\ell} + \beta_{\ell} - 1)(2n + \alpha_{\ell} + \beta_{\ell} - 2)},} \\
{\,\text{if}\, m = n - 1\,\text{and}\, n > 0,} \\
{\frac{1}{2} - \frac{\beta_{\ell}^{2} - \alpha_{\ell}^{2} - 2(\beta_{\ell} - \alpha_{\ell})}{2(2n + \alpha_{\ell} + \beta_{\ell})(2n + \alpha_{\ell} + \beta_{\ell} - 2)},} \\
{\,\text{if}\, m = n\,\text{and}\, n \geq 0,} \\
{\frac{(n + 1)(n + \alpha_{\ell} + \beta_{\ell} - 1)}{2(2n + \alpha_{\ell} + \beta_{\ell})(2n + \alpha_{\ell} + \beta_{\ell} - 1)},} \\
{\,\text{if}\, m = n + 1\,\text{and}\, n \geq 0,} \\
{0,} \\
{\text{otherwise},} \\
\end{array} \right.\ $$ which correspond to the coefficients of the three-term recurrence relation satisfied by the Jacobi Polynomials.

Matrix Expressions for the Transition Density
=============================================

It is computationally and notationally simpler to express the eigenfunctions of $\mathcal{L}_{\ell}$ and the transition density as products of matrices. In particular, we can express [equation (A.8)](#E53){ref-type="disp-formula"} as $$B_{\ell,n}(y) = e^{- {\overline{\sigma}}_{\ell}(y)/2}\mathbf{w}_{\ell,n}\mathbf{R}^{(\alpha_{\ell},\beta_{\ell})}(y),$$ where $$\mathbf{R}^{(\alpha,\beta)}(y) = {(R_{0}^{(\alpha,\beta)}(y),R_{1}^{(\alpha,\beta)}(y),\ldots)}^{T}$$ and we can express the vector $\mathbf{B}_{\ell}(y)$ of eigenfunctions as $$\mathbf{B}_{\ell}(y) = {(B_{\ell,0}(y),B_{\ell,1}(y),\ldots)}^{T}\, = e^{- {\overline{\sigma}}_{\ell}(y)/2}\mathbf{W}_{\ell}\mathbf{R}^{(\alpha_{\ell},\beta_{\ell})}(y),$$ where $$\mathbf{W}_{\ell} = {\lbrack\begin{array}{l}
\mathbf{w}_{\ell,0} \\
\mathbf{w}_{\ell,1} \\
 \vdots \\
\end{array}\rbrack}$$ is the matrix whose rows are the left eigenvectors of the matrix $\mathbf{M}_{\ell}$ in [equation (A.9)](#E54){ref-type="disp-formula"}.

Using [equations (A.5)](#E50){ref-type="disp-formula"} and (A.14), the transition density in a single epoch ℓ can then be expressed as the matrix product $$p_{\ell}(s,t;x,y) = \pi_{\ell}(y)\mathbf{B}_{\ell}^{T}(x)\mathbf{C}_{\ell}^{- 1}\mathbf{E}_{\ell}(t - s)\mathbf{B}_{\ell}(y),$$ where $$\mathbf{E}_{\ell}(t) = \text{diag}\left\{ e^{- \lambda_{\ell,0}t/2N_{\ell}},e^{- \lambda_{\ell,1}t/2N_{\ell}},\ldots \right\}$$ and $\mathbf{C}_{\ell} = \text{diag}\left\{ \left\langle B_{\ell,n},B_{\ell,n} \right\rangle_{\pi_{\ell}} \right\}_{n = 0}^{\infty}$. [@msw173-B48] showed that the matrix $\mathbf{C}_{\ell}$ in [equation (A.16)](#E61){ref-type="disp-formula"} can be expressed as $$\mathbf{C}_{\ell} = \mathbf{W}_{\ell}\mathbf{D}_{\ell}\mathbf{W}_{\ell}^{T},$$ where $$\mathbf{D}_{\ell} = \text{diag}\left\{ d_{0}^{(\alpha_{\ell},\beta_{\ell})},d_{1}^{(\alpha_{\ell},\beta_{\ell})},\ldots \right\}$$ and $$d_{i}^{(\alpha_{\ell},\beta_{\ell})} = \,\frac{\Gamma(i + \alpha_{\ell})\Gamma(i + \beta_{\ell})}{(2i + \alpha_{\ell} + \beta_{\ell} - 1)\Gamma(i + \alpha_{\ell} + \beta_{\ell} - 1)\Gamma(i + 1)}.$$

Thus, the transition density in a single epoch can be computed by constructing matrix $\mathbf{M}_{\ell}$, computing its eigenvectors $\mathbf{W}_{\ell}$ and eigenvalues $(\lambda_{\ell,0},\lambda_{\ell,1},\ldots)$, and plugging these into the components of [equation (A.16)](#E61){ref-type="disp-formula"}. In practice, because the matrix $\mathbf{M}_{\ell}$ has infinite dimension, we approximate it by truncating its dimensions at some large integer *M* yielding approximate eigenvectors $\left\{ {\widetilde{\mathbf{w}}}_{\ell,n} \right\}_{n = 0}^{M}$ and eigenvalues $\left\{ {\widetilde{\lambda}}_{\ell,n} \right\}_{n = 0}^{M}$. We also truncate the length of the vector $\mathbf{B}_{\ell}(y)$ at a large integer *N*. Although these truncations lead to approximate values of the transition density, the approximation can be made arbitrarily precise by taking $N \leq M$ to be sufficiently large.

Recursions for the Coefficients $\mathbf{a}_{k}$ and $\mathbf{b}_{k}$
=====================================================================

Discussion of the Problem
-------------------------

Here, we extend the HMM of [@msw173-B48] to accommodate populations of piecewise constant size. As we noted in the "Computing $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ under the Diffusion Approximation" section, the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ of the data under the diffusion model can be obtained using the equation $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = {\int_{y = 0}^{1}f_{K}}(y)dy,$$ where the quantity $f_{K}(y)$ is obtained by recursively evaluating [equations (16)](#E24){ref-type="disp-formula"} and [(17)](#E25){ref-type="disp-formula"}. Because $f_{k}(y)$ and $g_{k}(y)$ can be expressed as the series $f_{k}(y) = \pi_{\ell_{k}}(y)\mathbf{b}_{k}\mathbf{B}_{\ell_{k}}(y)$ and $g_{k}(y) = \pi_{\ell_{k}}(y)\mathbf{a}_{k}\mathbf{B}_{\ell_{k}}(y)$ ([eqs. 18](#E26){ref-type="disp-formula"} and [19](#E27){ref-type="disp-formula"}), determining $f_{k}(y)$ and $g_{k}(y)$ amounts to determining the coefficients $\mathbf{a}_{k}$ and $\mathbf{b}_{k}$. Thus, it is useful to develop analogs of the recursions [(15)](#E23){ref-type="disp-formula"} and [(16)](#E24){ref-type="disp-formula"} that apply to the coefficients themselves.

Equations for Propagating Coefficients
--------------------------------------

From [equation (16)](#E24){ref-type="disp-formula"}, it can be seen that obtaining $f_{k}(y)$ from $g_{k}(y)$ involves only multiplication by a polynomial in *y*. Thus, the formula for obtaining the coefficients $\mathbf{b}_{k}$ from the coefficients $\mathbf{a}_{k}$ does not depend on the population history and, therefore, it can be obtained from results in [@msw173-B48] who derived formulas for the recursion for the case of a population of constant size. However, the formula for obtaining $g_{k}(y)$ from $f_{k - 1}(y)$ ([eq. 17](#E25){ref-type="disp-formula"}) involves the transition probability $p_{\Theta}(t_{k - 1},t_{k};z,y)$, which depends on the population parameters Θ. Thus, it is necessary to account for the population history when computing the coefficients $\mathbf{a}_{k}$ from the coefficients $\mathbf{b}_{k - 1}$.

To obtain $\mathbf{a}_{k}$ from $\mathbf{b}_{k - 1}$, we first consider the more general problem of obtaining the generalized vector of coefficients $\mathbf{a}_{k}(t)$ from $\mathbf{b}_{k - 1}$, where $\mathbf{a}_{k}(t)$ is defined as the vector of coefficients of the expansion of the generalized density $g_{k}(y,t)$ defined by $$g_{k}(y,t)dy\,: = \mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:k - 1\rbrack} = o_{\lbrack 1:k - 1\rbrack},y \leq Y_{t} < y + dy \right\}\, = \pi_{\ell_{t}}(y)\mathbf{a}_{k}(t)\mathbf{B}_{\ell_{t}}(y),$$ i.e., the joint density of the observed data up to sample *k* -- 1 and the allele frequency at time *t*, where we assume $t_{k - 1} \leq t$ so that the time *t* at which $g_{k}(y,t)$ is evaluated occurs later than the time $t_{k - 1}$ at which $f_{k}(y)$ is evaluated. The density $g_{k}(y) = g_{k}(y,t_{k})$ defined in [equation (14)](#E22){ref-type="disp-formula"} is a special case of the generalized density $g_{k}(y,t)$ obtained when *t *=* t~k~*.

To obtain $\mathbf{a}_{k}(t)$ from $\mathbf{b}_{k - 1}$, there are two scenarios to consider: the case in which both $t_{k - 1}$ and *t* lie within the same epoch $\ell$ and the case in which $t_{k - 1}$ and *t* lie within distinct epochs. Our derivations of these separate cases provide the results necessary for step 2 of Procedure 2.

### The Case $\ell_{t_{k - 1}} = \ell_{t} = \ell$

If both $t_{k - 1}$ and *t* lie within the same epoch $\ell$, then the transition density is given by equation [(A.16)](#E61){ref-type="disp-formula"} and we have $$\pi_{\ell}(y)\mathbf{a}_{k}(t)\mathbf{B}_{\ell}(y)\, = g_{k}(y,t)\, = {\int_{0}^{1}f_{k - 1}}(z)p_{\ell}(t_{k - 1},t;z,y)dz\, = {\int_{0}^{1}\pi_{\ell}}(z)\mathbf{b}_{k - 1}\mathbf{B}_{\ell}(z)\pi_{\ell}(y)\mathbf{B}_{\ell}^{T}(z)\mathbf{C}_{\ell}^{- 1} \times \,\mathbf{E}_{\ell}(t - t_{k - 1})\mathbf{B}_{\ell}(y)dz\, = \pi_{\ell}(y)\mathbf{b}_{k - 1}{\lbrack{{\int_{0}^{1}\pi_{\ell}}(z)\mathbf{B}_{\ell}(z)\mathbf{B}_{\ell}^{T}(z)dz}\rbrack} \times \,\mathbf{C}_{\ell}^{- 1}\mathbf{E}_{\ell}(t - t_{k - 1})\mathbf{B}_{\ell}(y)\, = \pi_{\ell}(y)\mathbf{b}_{k - 1}\mathbf{E}_{\ell}(t - t_{k - 1})\mathbf{B}_{\ell}(y),$$ where the first equality follows from the definition of $g_{k}(y,t)$ (eq. B.2) and the second equality follows from [equation (17)](#E25){ref-type="disp-formula"}. In the fifth equality we have used the fact that ${\int_{0}^{1}\pi_{\ell}}(y)\mathbf{B}_{\ell}(y)\mathbf{B}_{\ell}^{T}(y)dy = \mathbf{C}_{\ell}$. Because the eigenfunctions $\left\{ B_{\ell,n}(y) \right\}_{n = 0}^{\infty}$ form a complete basis of the Hilbert space defined with respect to the inner product $\left\langle \cdot , \cdot \right\rangle_{\pi_{\ell}}$, the coefficients in the expansion on the left-hand side of [equation (B.3)](#E68){ref-type="disp-formula"} must equal those on the right-hand side. Thus, $$\mathbf{a}_{k}(t) = \mathbf{b}_{k - 1}\mathbf{E}_{\ell}(t - t_{k - 1}),\,\text{if}\,\ell_{t_{k - 1}} = \ell_{t}.$$

### The Case When $\ell_{t_{k - 1}} \neq \ell_{t}$

If the times $t_{k - 1}$ and *t* lie in different epochs, $\ell_{t_{k - 1}}$ and $\ell_{t}$, then the transition density is no longer given by [equation (A.16)](#E61){ref-type="disp-formula"}. Instead, we must use a formula for the transition density across multiple epochs of different sizes. [@msw173-B49] showed that if the allele frequency density $\rho_{\ell,s}(y)$ at time *s* in epoch $\ell$ is given by the expansion $$\rho_{\ell,s}(y) = \pi_{\ell}(y)\mathbf{r}_{\ell,s}\mathbf{B}_{\ell}(y),$$ where $\mathbf{r}_{\ell,s} = (r_{\ell,s,0},r_{\ell,s,1},\ldots)$ are the coefficients encoding the density at time *s* in the basis of the eigenfunctions $\left\{ \mathbf{B}_{\ell,n}(y) \right\}_{n = 0}^{\infty}$, then at time *t* in epoch $\ell + 1$, the allele frequency density is given by $\rho_{\ell + 1,t}(y) = \pi_{\ell + 1}(y)\mathbf{r}_{\ell + 1,t}\mathbf{B}_{\ell + 1}(y)$, where the coefficients $\mathbf{r}_{\ell + 1,t}$ are given by $$\mathbf{r}_{\ell + 1,t} = \mathbf{r}_{\ell,s}\mathbf{Z}_{\ell}(\tau_{\ell} - s;\zeta)\mathbf{E}_{\ell + 1}(t - \tau_{\ell}),$$ where $\tau_{\ell}$ is the time of the terminating boundary of epoch $\ell$, and $$\mathbf{Z}_{\ell}(\tau;\zeta)\, = \mathbf{E}_{\ell}(\tau)\mathbf{W}_{\ell}\mathbf{R}_{\ell}(\zeta)\mathbf{H}_{\ell,\ell + 1}(\zeta)\mathbf{R}_{\ell + 1}^{- 1}(\zeta)\mathbf{W}_{\ell + 1}^{- 1}.$$

In [equation (B.7)](#E72){ref-type="disp-formula"}, $\mathbf{R}_{\ell}(\zeta)$ and $\mathbf{H}_{\ell,\ell + 1}(\zeta)$ are given by $$\mathbf{R}_{\ell}(\zeta) = {\lbrack{\mathbf{R}^{(\alpha_{\ell},\beta_{\ell})}(\zeta_{0}),\mathbf{R}^{(\alpha_{\ell},\beta_{\ell})}(\zeta_{1}),\ldots}\rbrack},$$ where $\mathbf{R}^{\alpha,\beta}(y)$ is defined in [equation (A.13)](#E58){ref-type="disp-formula"} and $$\mathbf{H}_{\ell,\ell + 1}(\zeta)\, = \text{diag~}\left\{ {\frac{\pi_{\ell}(\zeta_{0})e^{- {\overline{\sigma}}_{\ell}(\zeta_{0})/2}}{\pi_{\ell + 1}(\zeta_{0})e^{- {\overline{\sigma}}_{\ell + 1}(\zeta_{0})/2}},} \right.\text{~~}\left. {\frac{\pi_{\ell}(\zeta_{1})e^{- {\overline{\sigma}}_{\ell}(\zeta_{1})/2}}{\pi_{\ell + 1}(\zeta_{1})e^{- {\overline{\sigma}}_{\ell + 1}(\zeta_{1})/2}},\ldots} \right\},$$ for an arbitrary collection of distinct values $\zeta = (\zeta_{0},\zeta_{1},\ldots) \in \lbrack 0,1\rbrack$. In practice, we take $\zeta$ to be the Chebyshev nodes ([@msw173-B49]).

By repeated application of [equation (B.6)](#E71){ref-type="disp-formula"}, it follows that if the coefficients $\mathbf{r}_{\ell_{s},s}$ encode the density $\rho_{s}(y)$ at time *s* in epoch $\ell_{s}$, then the coefficients $\mathbf{r}_{\ell_{t},t}$ encoding the density $\rho_{t}(y)$ at time *t* in epoch $\ell_{t} > \ell_{s}$ are given by $\mathbf{r}_{\ell_{t},t} = \mathbf{r}_{\ell_{s},s}\mathbf{F}(s,t;\zeta),$ where $$\mathbf{F}(s,t;\zeta)\, = \mathbf{Z}_{\ell_{s}}(\tau_{\ell_{s}} - s;\zeta)\left\lbrack {\prod\limits_{i = \ell_{s} + 1}^{\ell_{t} - 1}\mathbf{Z}_{i}(\tau_{i} - \tau_{i - 1};\zeta)} \right\rbrack \times \,\mathbf{E}_{\ell_{t}}(t - \tau_{\ell_{t} - 1}).$$ Moreover, if we define $\mathbf{r}_{\ell_{s},s}(x)$ to be the vector of coefficients encoding the density $\rho(y) = \delta(y - x)$, then it follows from [equation (B.10)](#E75){ref-type="disp-formula"} that the transition density $p_{\Theta}(s,t;x,y)$ for times *s *\<* t* lying in distinct epochs $\ell_{s} < \ell_{t}$ is given by $$\begin{array}{l}
{p_{\Theta}(s,t;x,y) = \pi_{\ell_{t}}(y)\mathbf{r}_{\ell_{s},s}(x)\mathbf{F}(s,t;\zeta)\mathbf{B}_{\ell_{t}}(y),} \\
{\,\text{if}\,\ell_{s} < \ell_{t}.} \\
\end{array}$$

For the initial condition $\rho_{\ell,s}(y) = \delta(y - x)$, it was shown in Proposition 1 of [@msw173-B48] that the coefficients $\mathbf{r}_{\ell_{s},s}(x)$ are given by $$\mathbf{r}_{\ell_{s},s}(x) = \left( {\frac{B_{\ell_{s},0}(x)}{\left\langle B_{\ell_{s},0},B_{\ell_{s},0} \right\rangle_{\pi_{\ell_{s}}}},\frac{B_{\ell_{s},1}(x)}{\left\langle B_{\ell_{s},1},B_{\ell_{s},1} \right\rangle_{\pi_{\ell_{s}}}},\ldots} \right)\, = \mathbf{B}_{\ell_{s}}{(x)}^{T}\mathbf{C}_{\ell_{s}}^{- 1},$$ yielding $$\begin{array}{l}
{p_{\Theta}(s,t;x,y) = \pi_{\ell_{t}}(y)\mathbf{B}_{\ell_{s}}{(x)}^{T}\mathbf{C}_{\ell_{s}}^{- 1}\mathbf{F}(s,t;\zeta)\mathbf{B}_{\ell_{t}}(y),\,} \\
{\,\text{if}\,\ell_{s} < \ell_{t},} \\
\end{array}$$ which is obtained by plugging [equation (B.12)](#E77){ref-type="disp-formula"} into [equation (B.11)](#E76){ref-type="disp-formula"}.

We can now plug [equation (B.13)](#E78){ref-type="disp-formula"} into [equation (17)](#E25){ref-type="disp-formula"} to obtain a relationship between $\mathbf{a}_{k}(t)$ and $\mathbf{b}_{k - 1}$ when times $t_{k - 1}$ and *t* lie in different epochs: $$\pi_{\ell_{t}}(y)\mathbf{a}_{k}(t)\mathbf{B}_{\ell_{t}}(y)\, = g_{k}(y,t)\, = {\int_{0}^{1}f_{k - 1}}(z)p_{\Theta}(t_{k - 1},t;z,y)dz\, = {\int_{0}^{1}\pi_{\ell_{k - 1}}}(z)\mathbf{b}_{k - 1}\mathbf{B}_{\ell_{k - 1}}(z)\pi_{\ell_{t}}(y)\mathbf{B}_{\ell_{k - 1}}{(z)}^{T}\,\mathbf{C}_{\ell_{k - 1}}^{- 1}\mathbf{F}(t_{k - 1},t;\zeta)\mathbf{B}_{\ell_{t}}(y)dz\, = \pi_{\ell_{t}}(y)\mathbf{b}_{k - 1}\ \left\lbrack {{\int_{0}^{1}\pi_{\ell_{k - 1}}}(z)\mathbf{B}_{\ell_{k - 1}}(z)} \right.\,\ \left. {\mathbf{B}_{\ell_{k - 1}}{(z)}^{T}dz} \right\rbrack\mathbf{C}_{\ell_{k - 1}}^{- 1}\mathbf{F}(t_{k - 1},t;\zeta)\mathbf{B}_{\ell_{t}}(y)\, = \pi_{\ell_{t}}(y)\mathbf{b}_{k - 1}\mathbf{F}(t_{k - 1},t;\zeta)\mathbf{B}_{\ell_{t}}(y),$$ where we have again used the fact that ${\int_{0}^{1}\pi_{\ell}}(z)\mathbf{B}_{\ell}(z)\mathbf{B}_{\ell}{(x)}^{T}dy = \mathbf{C}_{\ell}$. Finally, by the uniqueness of expansions in the Hilbert basis $\left\{ B_{\ell_{t},n} \right\}_{n = 0}^{\infty}$, we have $$\mathbf{a}_{k}(t) = \mathbf{b}_{k - 1}\mathbf{F}(t_{k - 1},t;\zeta),\,\text{if}\,\ell_{t_{k - 1}} \neq \ell_{t}.$$

The results derived in the "Equations for Propagating Coefficients" section provide the machinery necessary to propagate the coefficients $\mathbf{a}_{k}$ and $\mathbf{b}_{k}$ in the HMM over time. These results can now be used to compute the probability of observing a set of sampled allele frequencies under the diffusion model.

Derivation of Lemmas Necessary for Procedure 2
----------------------------------------------

We now obtain three lemmas that provide the steps in Procedure 2.Lemma B.3.1*If the initial frequency density* $\rho_{0}(\mathit{y})$ *at time* $\mathit{t}_{0} = 0$ *is* $\rho_{0}(\mathit{y}) = \delta(\mathit{y} - \mathit{x})$*, then the value of the initial vector* $\mathbf{b}_{0}$ *encoding the quantity* $\mathit{f}_{0}(\mathit{y})$ *is given by* $$\mathbf{b}_{0} = \left( {\frac{B_{\ell_{1},0}(x)}{c_{\ell_{1},0}},\frac{B_{\ell_{1},1}(x)}{c_{\ell_{1},1}},\ldots} \right) = \mathbf{C}_{\ell_{1}}^{- 1}\mathbf{B}_{\ell_{1}}(x),$$ where $\mathbf{B}_{\ell}(x)$ is given in equation (A.14) and $\mathbf{C}_{\ell}$ is the diagonal matrix given in [equation (A.18)](#E63){ref-type="disp-formula"}.ProofBecause $\mathbf{b}_{0}$ depends only on the parameters $\Theta_{\ell_{1}}$ in the first epoch, the proof of Lemma B.3.1 is the same whether we consider a population composed of a single epoch, or a population composed of multiple epochs. The equation for $f_{k}(y)$ ([eq. 18](#E26){ref-type="disp-formula"}) is the same as [equation 2.14](#E10){ref-type="disp-formula"} of [@msw173-B48]. Thus, the coefficients $\mathbf{b}_{k}$ in this paper correspond to the coefficients $\mathbf{b}_{k}$ in [@msw173-B48] who proved Lemma B.3.1 for the case of a population of constant size. Thus, the first equality in Lemma B.3.1 follows directly from Proposition 1 of [@msw173-B48]. The matrix representation in the second equality follows directly from the definitions of $\mathbf{C}_{\ell}$ and $\mathbf{B}_{\ell}(x)$. □Lemma B.3.2If the initial frequency density $\rho_{0}(y)$ at time $t_{0} = 0$ is unspecified and we instead assume a uniform prior on the allele frequency at the time of the first sampling event, then the coefficients $\mathbf{b}_{1}$ encoding the density at time t~1~ are given by $$\mathbf{b}_{1} = \frac{1}{\text{Beta}(o_{1} + 1,n_{1} - o_{1} + 1)B_{( - \sigma_{\ell_{1}},o_{1} + 1,n_{1} - o_{1} + 1),0}(0)} \times \,\mathbf{C}_{(\sigma_{\ell_{1}},o_{1} + 1,n_{1} - o_{1} + 1)}\mathbf{W}_{(\sigma_{\ell_{1}},o_{1} + 1,n_{1} - o_{1} + 1)} \times \,\mathbf{D}_{(\sigma_{\ell_{1}},o_{1} + 1,n_{1} - o_{1} + 1)}\mathbf{W}_{( - \sigma_{\ell_{1}},o_{1} + 1,n_{1} - o_{1} + 1)}^{T} \times \,\mathbf{Z}_{(\sigma_{\ell_{1}},o_{1} + 1,n_{1} - o_{1} + 1),(\sigma_{\ell_{1}},\alpha_{\ell_{1}},\beta_{\ell_{1}})}(0;\zeta),$$ where $\mathbf{C}_{(\sigma,\alpha,\beta)},\,\mathbf{D}_{(\sigma,\alpha,\beta)}$, and $\mathbf{W}_{(\sigma,\alpha,\beta)}$ are obtained by replacing $\sigma_{\ell},\,\alpha_{\ell}$, and $\beta_{\ell}$ in [equations (A.18)](#E63){ref-type="disp-formula"}, [(A.19)](#E64){ref-type="disp-formula"}, and [(A.15)](#E60){ref-type="disp-formula"} with σ, α, and β. Similarly, $\mathbf{Z}_{(\sigma,\alpha,\beta),(\widetilde{\sigma},\widetilde{\alpha},\widetilde{\beta})}(\tau;\zeta)$ is obtained by replacing $\sigma_{\ell},\,\sigma_{\ell + 1},\,\alpha_{\ell},\,\alpha_{\ell + 1},\,\beta_{\ell}$, and $\beta_{\ell + 1}$ in [equation (B.7)](#E72){ref-type="disp-formula"} with σ, $\widetilde{\sigma}$, α, $\widetilde{\alpha}$, β, and $\widetilde{\beta}$, respectively, $B_{(\sigma,\alpha,\beta)}(0)$ is the eigenfunction in [equation (B.23)](#E89){ref-type="disp-formula"} with $\sigma_{\ell},\,\alpha_{\ell}$, and $\beta_{\ell}$ replaced with σ, α, and β, and $\mathbf{W}_{(\sigma,\alpha,\beta)}$ is the matrix given in [equation (A.15)](#E60){ref-type="disp-formula"} whose rows are the left eigenvectors of $\mathbf{M}_{(\sigma,\alpha,\beta)}$, i.e., the matrix in [equation (A.9)](#E54){ref-type="disp-formula"} with $\sigma_{\ell},\,\alpha_{\ell}$, and $\beta_{\ell}$ replaced with σ, α, and β. Finally, $\text{Beta}(\alpha,\beta) = \Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha + \beta)$ is the normalizing constant of the beta distribution with parameters α and β.ProofIf the initial density $\rho_{0}(y)$ of the allele frequency is unspecified, then the density of $Y_{t_{1}}$ at the time *t*~1~ of the first sampling event, conditional on the observed number *o*~1~ of copies of allele *A* is given by Bayes Theorem: $$\left. f_{Y_{t_{1}}|O_{1}}(y \middle| o_{1}) = \frac{\mathbb{P}\left\{ O_{1} = o_{1} \middle| Y_{t_{1}} = y \right\} f_{Y_{1}}(y)}{\mathbb{P}\left\{ O_{1} = o_{1} \right\}}. \right.$$If we suppose that $f_{Y_{t_{1}}}(y)$ is uniform on the interval $\lbrack 0,1\rbrack$ then [equation (B.18)](#E83){ref-type="disp-formula"} becomes $$\left. f_{Y_{t_{1}}|O_{1}}(y \middle| o_{1})\, = \frac{\left( \begin{array}{l}
n_{1} \\
o_{1} \\
\end{array} \right)y^{o_{1}}{(1 - y)}^{n_{1} - o_{1}}}{{\int_{y = 0}^{1}\mathbb{P}}\left\{ O_{1} = o \middle| Y_{t_{1}} = y \right\} f_{Y_{t_{1}}}(y)dy}\, = \frac{\left( \begin{array}{l}
n_{1} \\
o_{1} \\
\end{array} \right)y^{o_{1}}{(1 - y)}^{n_{1} - o_{1}}}{{\int_{y = 0}^{1}\left( \begin{array}{l}
n_{1} \\
o_{1} \\
\end{array} \right)}y^{o_{1}}{(1 - y)}^{n_{1} - o_{1}}dy}\, = \frac{y^{o_{1}}{(1 - y)}^{n_{1} - o_{1}}}{{\int_{y = 0}^{1}y^{o_{1}}}{(1 - y)}^{n_{1} - o_{1}}dy}\, = \frac{y^{o_{1}}{(1 - y)}^{n_{1} - o_{1}}}{\text{Beta}(o_{1} + 1,n_{1} - o_{1} + 1)}, \right.$$ where, in the final equality, we have used the fact that $y^{o_{1}}{(1 - y)}^{n_{1} - o_{1}}$ is the unnormalized beta distribution with parameters $o_{1} + 1$ and $n_{1} - o_{1} + 1$. Thus, conditional on the first observation *o*~1~, we see that $\left. Y_{t_{1}} \middle| O_{1} \right.$ has a beta distribution with parameters $o_{1} + 1$ and $n_{1} - o_{1} + 1$.[@msw173-B48] showed in Appendix C of their paper that the distribution $\rho(y) = \frac{y^{\alpha - 1}{(1 - y)}^{\beta - 1}}{\text{Beta}(\alpha,\beta)}$ has coefficients given by $$\mathbf{b}_{\sigma,\alpha,\beta} = \frac{1}{\text{Beta}(\alpha,\beta)B_{( - \sigma,\alpha,\beta),0}(0)}\,\mathbf{C}_{(\sigma,\alpha,\beta)}\mathbf{W}_{(\sigma,\alpha,\beta)}\,\mathbf{D}_{(\sigma,\alpha,\beta)}\mathbf{W}_{( - \sigma,\alpha,\beta)}^{T}$$ in an epoch with selection parameter *σ* and mutation parameters *α* and *β*. Thus, using the change of basis in [equation (B.6)](#E71){ref-type="disp-formula"}, we see that the distribution $\rho(y) = \frac{y^{\alpha - 1}{(1 - y)}^{\beta - 1}}{\text{Beta}(\alpha,\beta)}$ has coefficients given by $\mathbf{b}_{\sigma,\alpha,\beta}\mathbf{Z}_{(\sigma,\alpha,\beta),(\sigma_{\ell},\alpha_{\ell},\beta_{\ell})}(0;\zeta)$ in the basis of eigenfunctions $\left\{ B_{\ell,n}(y) \right\}_{n = 0}^{\infty}$ corresponding to an epoch with parameters $\sigma_{\ell},\,\alpha_{\ell}$, and $\beta_{\ell}$, proving the result in [equation (B.17)](#E82){ref-type="disp-formula"}. □Lemma B.3.3Let $\mathbf{G}_{\ell},\,\mathbf{W}_{\ell},\,\mathbf{E}_{\ell}(t)$, and $\mathbf{F}(s,t;\zeta)$ denote the matrices defined in [equations (A.11)](#E56){ref-type="disp-formula"}, [(A.15)](#E60){ref-type="disp-formula"}, [(A.17)](#E62){ref-type="disp-formula"}, and [(B.10)](#E75){ref-type="disp-formula"}, respectively, where $\zeta = (\zeta_{0},\zeta_{1},\ldots)$ is a set of distinct values arbitrarily chosen such that $\left\{ \zeta_{0},\zeta_{1},\ldots \right\} \in \lbrack 0,1\rbrack$. Then the coefficient vectors $\mathbf{a}_{k}$ and $\mathbf{b}_{k}$ satisfy the recursive relationships $$\mathbf{b}_{k} = \mathbf{a}_{k}\mathbf{W}_{\ell_{k}}\mathbf{G}_{\ell_{k}}^{o_{k}}{(1 - \mathbf{G}_{\ell_{k}})}^{n_{k} - o_{k}}\mathbf{W}_{\ell_{k}}^{- 1},\,$$ $$a_{k} = \begin{cases}
{\mathbf{b}_{k - 1}\mathbf{E}_{l_{k}}(t_{k} - t_{k - 1})} & {if\,\ell_{k - 1} = \ell_{k},} \\
{\mathbf{b}_{k - 1}\mathbf{F}(t_{k - 1},t_{k};\zeta)} & {otherwise,} \\
\end{cases}$$ where $\mathbf{W}_{\ell}^{- 1} = \mathbf{D}_{\ell}\mathbf{W}_{\ell}^{T}\mathbf{C}_{\ell}^{- 1}$.ProofThe relationship in [equation (B.21)](#E87){ref-type="disp-formula"} is obtained immediately by setting *t* = *t*~*k*~ in [equations (B.4)](#E69){ref-type="disp-formula"} and [(B.15)](#E80){ref-type="disp-formula"}, which follows because $\mathbf{a}_{k}(t_{k}) = \mathbf{a}_{k}$. The relationship in [equation (B.20)](#E86){ref-type="disp-formula"} does not depend on the population parameters Θ; therefore, [equation (B.20)](#E86){ref-type="disp-formula"} is the same as that derived in [@msw173-B48], who considered a population of constant size (see [@msw173-B48], Theorem 2). □Lemma B.3.4The probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}$ of observing the allele counts $o_{\lbrack 1:K\rbrack}$, given the population parameters Θ is $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\} = \frac{c_{\ell_{K},0}}{B_{\ell_{K},0}(0)}b_{K,0},$$ where $c_{\ell,0} = {\lbrack\mathbf{C}_{\ell}\rbrack}_{0,0}$ is element 0, 0 of the matrix $\mathbf{C}_{\ell}$ in [equation (A.18)](#E63){ref-type="disp-formula"} and $$B_{\ell,0}(0) = \sum\limits_{m = 0}^{\infty}{( - 1)}^{m}{\lbrack\mathbf{W}_{\ell}\rbrack}_{0,m}\frac{\Gamma(m + \alpha_{\ell})}{\Gamma(m + 1)\Gamma(\alpha_{\ell})}.$$The quantity ${\lbrack\mathbf{W}_{\ell}\rbrack}_{i,j}$ in [equation (B.23)](#E89){ref-type="disp-formula"} is element i, j of the matrix $\mathbf{W}_{\ell}$ given in [equation (A.15)](#E60){ref-type="disp-formula"}.Proof[Equation (B.22)](#E88){ref-type="disp-formula"} can be obtained by integrating over the joint density $f_{K}(y)$ of the data $O_{\lbrack 1:K\rbrack}$ and the allele frequency $Y_{t_{K}}$ at the final sampling time: $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}\left. \, = {\int_{0}^{1}f_{K}}(y)dy\, = {\int_{0}^{1}{\sum\limits_{n = 0}^{\infty}b_{K,n}}}\pi_{\ell_{K}}(y)B_{\ell_{K},n}(y)dy\, = \sum\limits_{n = 0}^{\infty}b_{K,n}{\int_{0}^{1}\pi_{\ell_{K}}}(y)B_{\ell_{K},n}(y)dy\, = \sum\limits_{n = 0}^{\infty}b_{K,n}{\int_{0}^{1}\pi_{\ell_{K}}}(y)B_{\ell_{K},n}(y \right)\frac{B_{\ell_{K},0}(y)}{B_{\ell_{K},0}(0)}dy\, = b_{K,0}\frac{c_{\ell_{K},0}}{B_{\ell_{K},0}(0)},$$ where $c_{\ell_{K},0} = {\lbrack\mathbf{C}_{\ell_{K}}\rbrack}_{0,0} \equiv \left\langle B_{\ell_{K},0},B_{\ell_{K},0} \right\rangle_{\pi_{\ell_{K}}}$. In the fourth equality we have used the fact that $B_{\ell,0}(y) = B_{\ell,0}(0)$ is a constant function in *y*. To see why $B_{\ell,0}(y)$ is constant, note that the eigenvalues $\lambda_{\ell,0},\lambda_{\ell,1},\ldots$ are non-negative and strictly increasing. Thus, all terms in [equation (A.5)](#E50){ref-type="disp-formula"} must vanish in the limit $\left. s\rightarrow - \infty \right.$, except possibly the term *n* = 0. Because $p_{\ell}(s,t;x,y)$ approaches the stationary density in the limit $\left. s\rightarrow - \infty \right.$, it must be the case that $\lambda_{\ell,0} = 0$, so at least one term does not vanish. Thus, we have $$\left. \underset{s\rightarrow - \infty}{\lim}p_{\theta}(s,t;x,y) = \pi_{\ell}(y \right)\frac{B_{\ell,0}(x)\, B_{\ell,0}(y)}{\left\langle B_{\ell,0},B_{\ell,0} \right\rangle_{\pi_{\ell}}} \propto \pi_{\ell}(y),$$ where we have used the fact that $\pi_{\ell}(y)$ is proportional to the stationary density of the diffusion equation in Epoch $\ell$. It follows from [equation (B.25)](#E91){ref-type="disp-formula"} that $B_{\ell,0}(y)$ is constant. Thus, we obtain the result, proving [equation (B.22)](#E88){ref-type="disp-formula"}. [Equation (B.23)](#E89){ref-type="disp-formula"} follows directly from the proof of Proposition 3 in [@msw173-B48]. □

Conditional Probabilities: The Wright--Fisher Model
===================================================

Under the Wright--Fisher model, the probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}$ of the observed allele counts, conditional on the event *S~K~* that allele *A* is segregating in the final sample can be computed using the fact that $$\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}\, = \sum\limits_{j = 0}^{2N_{t_{K}}}\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack},c_{t_{K}} = j \middle| S_{K} \right\}\, = \sum\limits_{j = 0}^{2N_{t_{K}}}\frac{\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \middle| O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack},c_{t_{K}} = j \right\}}{\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \right\}} \times \,\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack},c_{t_{K}} = j \right\}\, = \frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \right\}} \times \,\sum\limits_{j = 0}^{2N_{\ell_{t_{K}}}}\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack},c_{t_{K}} = j \right\}\, = \frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \right\}}\sum\limits_{i = 0}^{2N_{t_{K}}}\mathbf{v}_{K,i},$$ where the third equality in [equation (C.1)](#E92){ref-type="disp-formula"} follows from the fact that the conditional probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \middle| O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack},c_{t_{K}} = j \right\}$ depends only on the allele count *o~K~* and the final equality in [equation (C.1)](#E92){ref-type="disp-formula"} follows from the definition of $\mathbf{v}_{k}$. The probability $\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\}$ in [equation (C.1)](#E92){ref-type="disp-formula"} is given by $$\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\} = \begin{cases}
{1,} & {\,\text{if}\, 1 \leq o_{K} < n_{K},} \\
{0,} & {\,\text{otherwise}\,} \\
\end{cases}\ $$ and the probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \right\}$ is given by $$\mathbb{P}_{\Theta,\mathcal{W}}\left\{ S_{K} \right\}\, = \sum\limits_{i = 0}^{2N_{t_{K}}}\mathbb{P}\left\{ S_{K} \middle| C_{t_{K}} = i \right\}\mathbb{P}_{\Theta,\mathcal{W}}\left\{ C_{t_{K}} = i \right\}\, = \sum\limits_{i = 0}^{2N_{t_{K}}}{\ \left\lbrack {1 - \mathbb{P}\left\{ O_{K} = 0 \middle| C_{t_{K}} = i \right\} - \ } \right.}\mathbb{P}\left\{ O_{K} = n_{K} \middle| C_{t_{K}} = i \right\}\rbrack \times \,\mathbb{P}_{\Theta,\mathcal{W}}\left\{ C_{t_{K}} = i \right\}\, = \sum\limits_{i = 0}^{2N_{t_{K}}}\left\lbrack {1 - \left( {1 - \frac{i}{2N_{t_{K}}}} \right)^{n_{K}} - \left( \frac{i}{2N_{t_{K}}} \right)^{n_{K}}} \right\rbrack \times \,\mathbb{P}_{\Theta,\mathcal{W}}\left\{ c_{t_{K}} = i \right\},$$ where, as before, $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ C_{t_{K}} = i \right\}$ is given by the *i*th element of $\mathbf{d}_{t} = \mathbf{d}_{0}\prod_{t = 1}^{t_{K}}T_{t - 1,t}$.

Note that it is easy to condition on other configurations of the final sample using a procedure similar to that used to derive [equation (C.1)](#E92){ref-type="disp-formula"}. For example, for the event F~K~ that allele A is segregating or fixed in the final sample, which we consider in the "Simulations Conditioning on Segregation or Fixation" section, the probability $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| F_{K} \right\}$ is given by $$\mathbb{P}_{\Theta,\mathcal{W}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| F_{K} \right\}\, = \frac{\mathbb{P}\left\{ F_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{W}}\left\{ F_{K} \right\}}\sum\limits_{i = 0}^{2N_{t_{K}}}\mathbf{v}_{K,i},$$ where $$\begin{array}{l}
{\mathbb{P}\left\{ F_{K} \middle| O_{K} = o_{K} \right\}} \\
{= \begin{cases}
{1,} & {\,\text{if}\, 1 \leq o_{K} \leq n_{K},} \\
{0,} & \text{otherwise} \\
\end{cases}} \\
\end{array}$$ and $$\mathbb{P}_{\Theta,\mathcal{W}}\left\{ F_{K} \right\}\, = \sum\limits_{i = 0}^{2N_{t_{K}}}\left\lbrack {1 - \mathbb{P}\left\{ O_{K} = 0 \middle| C_{t_{K}} = i \right\}} \right\rbrack\mathbb{P}_{\Theta,\mathcal{W}}\left\{ C_{t_{K}} = i \right\}\, = \sum\limits_{i = 0}^{2N_{t_{K}}}\left\lbrack {1 - \left( {1 - \frac{i}{2N_{t_{K}}}} \right)^{n_{K}}} \right\rbrack\mathbb{P}_{\Theta,\mathcal{W}}\left\{ c_{t_{K}} = i \right\}.$$

Other probabilities can be obtained in a similar fashion.

Conditional Probabilities: Diffusion Model
==========================================

Under the diffusion approximation, the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}$ of the observed allele counts conditional on the event *S~K~* that allele *A* is segregating in the final sample can be computed using the fact that $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}\, = {\int_{y = 0}^{1}\mathbb{P}_{\Theta,\mathcal{D}}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack},Y_{t_{K}} = y \middle| S_{K} \right\} dy\, = {\int_{y = 0}^{1}\frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{D}}\left\{ S_{K} \right\}}}f_{K}(y)dy\, = \frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{D}}\left\{ S_{K} \right\}}{\int_{0}^{1}f_{K}}(y)dy\, = \frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{D}}\left\{ S_{K} \right\}}\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \right\}\, = \frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\}}{\mathbb{P}_{\Theta,\mathcal{D}}\left\{ S_{K} \right\}}\frac{c_{\ell_{K},0}}{B_{\ell_{K},0}(0)}b_{K,0},$$ where the second equality follows from the fact that the conditional probability $\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K},Y_{t_{K}} = y \right\}$ depends only on the allele count *o~K~* in the final sample, and the final equality follows from [equation (B.22)](#E88){ref-type="disp-formula"}.

The probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ S_{K} \right\}$ can be computed as $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ S_{K} \right\}\, = 1 - \mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{K} = 0 \right\} - \mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{K} = n_{K} \right\}.$$

In [equation (D.2)](#E99){ref-type="disp-formula"}, the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{K} = j \right\}$ can be found easily by noting that if the only sampling time is *t~K~*, at which *O~K ~*=*~ ~j* lineages are observed, then the probability computed using Procedure 2 is precisely the probability $\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{K} = j \right\}$.

Consider the problem in which the only sampling occurs at time *t~K~* and denote the coefficient vectors for this related problem by ${\widetilde{\mathbf{a}}}_{k}$ and ${\widetilde{\mathbf{b}}}_{k}$. Then, by [equation (B.22)](#E88){ref-type="disp-formula"}, we see that $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{K} = j \right\} = \frac{c_{\ell_{K},0}}{B_{\ell_{K},0}(0)}{\widetilde{b}}_{K,0}(j),$$ where ${\widetilde{b}}_{K}(j)$ is obtained by computing the steps in Procedure 2. In Step 1, we compute $${\widetilde{\mathbf{b}}}_{0} = \mathbf{b}_{0},$$ which follows because the initial vector $\mathbf{b}_{0}$ depends only on the initial frequency. In Step 2, we compute $${\widetilde{\mathbf{a}}}_{K} = \ \begin{cases}
{{\widetilde{\mathbf{b}}}_{0}\mathbf{E}_{\ell_{1}}(t_{K}),} & {\,\text{if}\,\ell_{t_{K}} = 1,} \\
{{\widetilde{\mathbf{b}}}_{0}\mathbf{F}(0,t_{K};\zeta),} & {\text{otherwise},} \\
\end{cases}$$ which follows because the coefficients are propagated directly from time $t_{0} = 0$ to time *t~K~*. Finally, in Step 3 we have $${\widetilde{b}}_{K}(j) = {\widetilde{\mathbf{a}}}_{K}\mathbf{W}_{\ell_{K}}\mathbf{G}_{\ell_{K}}^{j}{(1 - \mathbf{G}_{\ell_{K}})}^{n_{K} - j}\mathbf{W}_{\ell_{K}}^{- 1}.$$

Combined together, [equations (D.4)](#E101){ref-type="disp-formula"}, [(D.5)](#E102){ref-type="disp-formula"}, and [(D.6)](#E103){ref-type="disp-formula"} yield $${\widetilde{\mathbf{b}}}_{K}(j)\  = \ \left\{ \begin{array}{l}
{\mathbf{b}_{0}\mathbf{E}_{\ell_{1}}(t_{K})\mathbf{W}_{\ell_{K}}\mathbf{G}_{\ell_{K}}^{j}{(1 - \mathbf{G}_{\ell_{K}})}^{n_{K} - j}\mathbf{W}_{\ell_{K}}^{- 1},} \\
{\,\text{if}\,\ell_{t_{K}} = 1,} \\
{\mathbf{b}_{0}\mathbf{F}(0,t_{K};\zeta)\mathbf{W}_{\ell_{K}}\mathbf{G}_{\ell_{K}}^{j}{(1 - \mathbf{G}_{\ell_{K}})}^{n_{K} - j}\mathbf{W}_{\ell_{K}}^{- 1},} \\
{\text{otherwise}.} \\
\end{array} \right.$$

Plugging [equations (D.2)](#E99){ref-type="disp-formula"} and [(D.3)](#E100){ref-type="disp-formula"} into [equation (D.1)](#E98){ref-type="disp-formula"} gives $$\mathbb{P}_{\Theta,\mathcal{D}}\left\{ O_{\lbrack 1:K\rbrack} = o_{\lbrack 1:K\rbrack} \middle| S_{K} \right\}\, = \frac{\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\} c_{\ell_{K},0}b_{K,0}}{B_{\ell_{K},0}(0) - c_{\ell_{K},0}{\widetilde{b}}_{K,0}(0) - c_{\ell_{K},0}{\widetilde{b}}_{K,0}(n_{K})},$$ where $$\mathbb{P}\left\{ S_{K} \middle| O_{K} = o_{K} \right\} = \begin{cases}
{1,} & {\,\text{if}\, 1 \leq o_{K} < n_{K},} \\
{0,} & {\text{otherwise}.} \\
\end{cases}$$

Note that it is easy to condition on other configurations of the final sample by computing the probabilities $\mathbb{P}\left\{ V_{K} \middle| O_{K} = o_{k} \right\}$ and $\mathbb{P}_{\Theta,\mathcal{W}}\left\{ V_{K} \right\}$ for some other event *V~K~*.

[^1]: **Associate editor:** Ryan Hernandez
