EWF: simulating exact paths of the Wright–Fisher diffusion

Abstract Motivation The Wright–Fisher diffusion is important in population genetics in modelling the evolution of allele frequencies over time subject to the influence of biological phenomena such as selection, mutation and genetic drift. Simulating the paths of the process is challenging due to the form of the transition density. We present EWF, a robust and efficient sampler which returns exact draws for the diffusion and diffusion bridge processes, accounting for general models of selection including those with frequency dependence. Results Given a configuration of selection, mutation and endpoints, EWF returns draws at the requested sampling times from the law of the corresponding Wright–Fisher process. Output was validated by comparison to approximations of the transition density via the Kolmogorov–Smirnov test and QQ plots. Availability and implementation All softwares are available at https://github.com/JaroSant/EWF. Supplementary information Supplementary data are available at Bioinformatics online.


Introduction
The Wright-Fisher diffusion is a central model for the temporal fluctuation of allele frequencies in a large population evolving under random mating and in the presence of mutation and selection.Despite its importance, it remains difficult to work with from a computational perspective, both in the absence of selection (where the transition density admits an infinite series expansion) and the non-neutral case (where the corresponding infinite series expansion has intractable terms).Additionally, in a diallelic model the diffusion lives on the bounded interval [0, 1] and thus even simple approximate sampling techniques such as the Euler-Maruyama scheme require sophisticated modifications to respect its boundary behaviour (Dangerfield et al., 2012).
Existing approaches in the literature have tackled this by resorting to a combination of discretisation and numerical approximation, e.g.solving the Kolmogorov backwards equation numerically (Bollback et al., 2008;Malaspinas et al., 2012), approximating through more tractable processes (Mathieson and McVean, 2013), truncating a spectral expansion of the transition density (Steinrücken et al., 2016), and using Riemann sum approximations (Schraiber et al., 2016), all of which induce a bias which is hard to quantify.
In some cases, exact sampling routines making use of rejection sampling are available.This class of techniques has been extended to certain variants of the Wright-Fisher diffusion: Jenkins and Spanò (2017) showed that neutral Wright-Fisher diffusion paths and bridges can be simulated exactly via simulation techniques tailored for infinite series, and that neutral paths are the natural proposal mechanism for simulating non-neutral paths by rejection.Their work assumes that the mutation parameters are strictly positive and the endpoints for both the diffusion and diffusion bridge lie in the interior of [0,1].The case of diffusion bridges that start and end at 0 was tackled by Griffiths et al. (2018), but several other combinations of startpoint, endpoint, and parameters remain unaddressed.Moreover, no simulation package implementing all of the cases of interest exists.
We present EWF, a C++ package producing exact draws from both neutral and non-neutral Wright-Fisher diffusions.The method properly accounts for all types of boundary (entrance, reflecting, and absorbing), incorporates a wide class of selection models, and allows for arbitrary endpoints, substantially extending previous work by Jenkins and Spanò (2017); Griffiths et al. (2018).These new theoretical details can be found in the accompanying supplement.Additionally, EWF preserves accuracy over long times, in contrast to Euler-Maruyama type schemes where errors accumulate over the simulated path.

Models
Consider the two-allele non-neutral Wright-Fisher diffusion (X t ) t≥0 with mutation parameter θ = (θ 1 , θ 2 ), which is given by the solution to the following stochastic differential equation for t ≥ 0 with X 0 ∈ [0, 1], and η(x) = n i=0 a i x i for n finite (e.g. for genic selection η(x) = 1, and for diploid selection η(x) = h + x(1 − 2h) with h the dominance parameter).When the mutation parameter θ has positive entries, the corresponding neutral (i.e.σ = 0) transition density can  2009) simulated using EWF (note that the observed frequencies (black crosses) are assumed to be exact observations of the underlying diffusion).Simulations used the inferred selection coefficient s = 0.0007 with a consensus effective population size N e = 10, 000 (Ludwig et al., 2009;Malaspinas et al., 2012;Schraiber et al., 2016), giving σ = 2N e s = 14.We used θ = 0 and a generation time of 5 years.
The implementation was tested extensively and validated through a combination of QQ plots and the Kolmogorov-Smirnov test (see Supplementary Information Section 7).An example is shown in Fig. 1.

Discussion
EWF provides a robust, efficient, and exact sampling routine to target a wide family of Wright-Fisher diffusions featuring a broad class of selective regimes, any mutation parameters, and any start/end points.The implementation can be used as a stand-alone package, or incorporated into simulation-based inference pipelines from time series allele frequency data.This is particularly useful in view of the recent increase in availability of such data (Wutke et al., 2016;Fages et al., 2019).
and θ 2 + m − l, and , with Γ(•) denoting the gamma function.We point out that {q θ m (t)} m∈N correspond to the transition probabilities of the number of lineages in Kingman's coalescent (which is the moment dual to the Wright-Fisher diffusion), such that q θ m (t) is the probability that m lineages survive up to time t when one starts with an infinite number of lineages at time 0. For more details, we refer the interested reader to Griffiths (1979); Tavaré (1984).The inclusion of the mutation parameters on the LHS of (2) makes explicit the dependence of the transition density on these quantities, however in an effort to reduce on encumbrance, we shall suppress this notation henceforth and simply write p(x, y; t) for the transition density of the diffusion, with the specific mutation regime being considered specified exogenously.

Neutral diffusion with one sided mutation
For θ = (0, θ), the diffusion is absorbed upon hitting 0 and the transition density can be expressed as where δ 0 (y) denotes a point mass at 0 and represents the case when the diffusion is absorbed at 0. In cases like this we reinterpret 'density' appropriately, with respect to a dominating measure containing both a Lebesgue component and an atom at each of 0 and 1.
If we condition on the event {t < T 0 }, standard conditional probability gives us that the transition density of the diffusion conditioned on non-absorption until time t is given by p(x, y; t) = p(x, y; t) for y ∈ (0, 1], where we use the notation p(•, •; •) to make explicit the fact that this is the transition density of the conditioned diffusion process.Additionally, we have that and we note that the contributions from m = 0 above are missing as the corresponding beta density collapses to a point mass at 0. Thus for x, y ∈ (0, 1] we have (5) For small x, we have the following leading order expansion in x p(x, y; and note further ( 4) is also of leading order x for x small.Thus upon taking the limit x → 0 in (5) we get that Putting all of the above together we get that the conditioned diffusion has transition density given by p(x, We point out that as the diffusion is conditioned on avoiding 0, there will always be at least one surviving lineage in the moment-dual Kingman coalescent, and thus the index for m starts at 1.

Diffusion without mutation
If θ = 0, then the diffusion is absorbed upon hitting either boundary, and the corresponding transition density is given by Conditioning the diffusion on remaining inside the interior of [0, 1], and again employing a leading order analysis of the resulting numerator and denominator allows us to conclude that the transition density in this case is given by p Note that as θ = 0 and we are conditioning on non-absorption, the indices m and d are now forced to start from 2. This follows from the fact that the derivations performed above assume the starting point x to be within (0, 1) and subsequently send x to the corresponding boundary from within the interior of (0, 1), which corresponds to starting the diffusion arbitrarily close to the boundary.Thus at all times there is a fraction x of the population having one type, with the other fraction 1 − x having the other, neither of which can be lost by mutation.

Transition densities for neutral Wright-Fisher diffusion bridges
We now derive the density of a point y ∈ [0, 1] sampled at time s ∈ (0, t) from the law of a Wright-Fisher diffusion bridge started at x at time 0 and ending at z at time t.In addition to considering each mutation regime separately, we further split our considerations based on the values the start and end points x and z assume.As in the diffusion case, we derive the relevant expressions in the case θ = (0, θ), as the other cases (θ = (0, 0) or θ = (θ, 0)) follow using similar arguments.We further consider both cases when (i) the bridge is allowed to be absorbed at any time point within the time interval (0, t), and (ii) the bridge is conditionally non-absorbing: X s ∈ (0, 1) for all s ∈ (0, t).We make use of the following short-hand notation for the different possible end-point combinations.
We further introduce a letter at the front of each of the above to differentiate between the cases θ = 0 ('Z' for zero), θ = (0, θ) ('O' for one sided), and θ with strictly positive entries ('P' for strictly positive).
Before proceeding with deriving the transition densities for all the above outlined cases, observe that the transition density for a Wright-Fisher diffusion bridge started from x ∈ [0, 1] at time 0, ending at z ∈ [0, 1] at time t and sampled at time s can be factorised as follows for y ∈ [0, 1]: p x,z;t (y; s) = p(x, y; s)p(y, z; t − s) p(x, z; t) , where again, for simplicity the dependence of (11) on the mutation parameters is omitted from the notation.
2.1 Neutral diffusion bridge with one sided mutation θ = (0, θ) We start by noting that if the diffusion bridge is allowed to be absorbed at 0 at any time within the interval (0, t), then the only cases of interest are when the left endpoint x ∈ (0, 1], for otherwise the bridge stays at 0. Additionally if z ∈ (0, 1], the bridge could not have been absorbed within the time interval (0, t), and is therefore equivalent to conditioning it on nonabsorption (which shall be tackled shortly).Thus we take x ∈ (0, 1) and z = 0, substitute (3) into (11), and re-group terms to get that where B(•, •) is the beta function.
To derive the transition density when z ∈ (0, 1], we first point out that conditioning a diffusion (or conditioning a diffusion bridge) on non-absorption is a special case of taking an h-transform for said process (see for instance and thus the distribution of a diffusion bridge conditioned on non-absorption is the same as that of the corresponding unconditioned process.We therefore need not differentiate between the transition density of the conditioned or unconditioned diffusion bridge, and simply use p x,z;t (y; s) throughout.
Expanding (11) for x, z ∈ (0, 1] gives When x = 0, we make use of (6) in both the numerator and denominator above, and subsequently take the limit as x → 0, to arrive at The above expression can further be used to derive the expression when z = 0 by taking leading order terms in z and taking the limit z → 0, giving As previously mentioned, the case θ = (θ, 0) follows from the above by considering the symmetric map x → 1 − x.

Neutral diffusion bridge with no mutation
We can replicate all of the above arguments for when θ = 0 to get that if x ∈ (0, 1) and z = 0, for y ∈ [0, 1) we have whilst if x ∈ (0, 1) and z = 1, we get for y ∈ (0, 1] Note that if z = 0, then we cannot have y = 1 and similarly if z = 1, y cannot be equal to 0.
Computing the transition densities conditioned on non-absorption can be done as in the onesided mutation case illustrated above, by following the same arguments.
The resulting expressions for the conditioned diffusion bridges under all three mutation regimes can be found below (recall the notation in Table 2).

Sampling schemes
We now detail how to obtain sample draws from the above transition densities for both the diffusion and diffusion bridge case.
In this case, Algorithm 1 in Jenkins and Spanò (2017) can be easily adapted to sample from (3): The only modification to Algorithm 1 in Jenkins and Spanò (2017) is the sampling procedure in step 3, where the outcome L = 0 encodes the event when the diffusion is absorbed before the sampling time.A similar strategy allows for draws from (9), where additionally if L = m, then in step 3 we return Y = 1.
For the case when the diffusion is conditioned on non-absorption, both expressions on the RHS of ( 8) are mixtures of beta distributions, with the weights forming a probability mass function (pmf) on N. When the starting point x is set to 0, one can return a draw Y from the law of the corresponding diffusion process sampled at time t, by Step 2 is straightforward, whilst for step 1 the 'alternating series trick' can be employed-this technique requires access to a pair of monotonic sequences of upper and lower bounds for terms in the numerator and denominator, both converging to their exact values.This is immediate for the numerator (Proposition 1 in Jenkins and Spanò (2017)), whilst for the denominator we modify slightly the arguments present in Proposition 3 in Jenkins and Spanò (2017) (see Section 5 for further details).
A similar sampling scheme can be used for drawing samples from the law of the diffusion started from x ∈ (0, 1], where once again appropriate monotonic upper and lower bounds can be constructed for both numerator and denominator (see Section 5 for more details).
The above can be replicated and suitably tweaked to return samples from (10), where an additional scheme is needed to deal with the case x = 1.

Sampling from the law of the diffusion bridge
Once again we start by considering the case when the bridge is allowed to be absorbed at the boundary within the time interval (0, t), and the mutation parameter is given by θ = (0, θ).
To sample from ( 12), we follow an approach similar to that illustrated above for the diffusion case.Recall that we need only focus on the case when z = 0, for otherwise the bridge cannot have been absorbed during the time interval (0, t) and thus is equivalent to conditioning on non-absorption (for which an appropriate sampling scheme will be provided shortly).Observe that the RHS of ( 12) can be viewed as a mixture of beta distributions, with the mixture weights defining a pmf on a subspace of N 3 (for more details, refer to (Jenkins and Spanò, 2017, Section 3)).Individual monotonic upper and lower bounds can be constructed for {q θ m (s)} m∈N , {q θ k (t − s)} k∈N and ∞ d=1 q θ d (t)(1 − x) d (see Section 5 for full details with regards to this last quantity), and subsequently these can be put together to obtain monotonic upper and lower bounds on the {p m,k,l } m,k,l∈N .Thus the alternating series trick lends itself to return a draw (M, K, L) ∼ {p m,k,l } m,k,l,∈N , and we use this to draw the relevant sample diffusion bridge point: A similar scheme can be derived for the case θ = (θ, 0) by symmetric arguments, whereas for θ = 0 the above can be replicated with the only significant difference being that if L = m, then the routine returns Y = 1 in step 2.
We now turn to the case when the diffusion bridge is conditioned on not being absorbed within the time interval (0, t).Corollary 2 in Griffiths et al. (2018) gives us that Wright-Fisher diffusion bridges with mutation parameters either θ = 0 or θ = (0, θ) are equal (in distribution) to Wright-Fisher bridges with mutation parameters θ = (2, 2) or θ = (2, θ) respectively.Thus from now on we shall focus our attention solely on the case when θ 1 , θ 2 > 0.
The strategy will be very close to the one developed above and based on the method found in (Jenkins and Spanò, 2017, Section 3).As in the unconditioned bridge case, the diffusion bridge densities (PA1)-(PC3) can be viewed as mixtures of beta distributions, where the mixture weights now define a pmf on subspaces of N 4 and whose exact form depends on the particular density being considered.
As diffusion bridges are invariant under time reversal, a diffusion bridge that goes from x to y in time s and then proceeds to terminate at z at time t has the same law as a diffusion bridge that starts at z, proceeds to y at time t − s and ends at x at time t.This, coupled with symmetric arguments allows us to sample from the various transition densities (PA1)-(PC3) using just four different schemes, which we group as follows: 1. Start and endpoints are the same (i.e.equations (PA1) and (PB3)).
3. z is in the interior of [0, 1], and the starting point is at one of the boundary points (i.e.equations (PA2), (PB2), (PC1) and (PC3) -note that for (PC1) and (PC3) we make use of time reversal).
Using the above groupings, it remains to show that the resulting four different transition densities consist of mixture weights {p m,k,l,j } m,k,l,j∈N for which one can obtain monotonic sequences of upper and lower bounds.Again constructing these quantities for the numerator is straightforward, whereas the denominator is tackled in Section 5 (by suitably modifying Proposition 4 from Jenkins and Spanò (2017)).
We point out that for both the diffusion and diffusion bridge case, numerical instabilities present when computing contributions to the infinite series representation of the probabilities {q θ m (t)} m∈N for small time increments prompt the use of approximations for these quantities.For more details, please refer to Section 6.

Simulation of non-neutral paths
As observed in (Jenkins and Spanò, 2017, Section 5), the neutral Wright-Fisher process can be used as a proposal distribution in an appropriate rejection sampler to returns exact draws from a non-neutral process.We give a brief overview for completeness.
Denote by WF x 0 σ,θ the law induced by the solution X T := (X t ) T t=0 to the SDE given by equation (1) in the main paper on the space of continuous functions mapping [0, T ] into [0, 1] for some fixed time T , and by WF x 0 0,θ the corresponding neutral law.The Radon-Nikodym derivative between these two laws is given by where Ã(x) := (σ/2) x 0 η(z)dz with Ã(x) ≤ Ã+ for any x ∈ [0, 1], and Observe that ϕ(x) is a polynomial in x (in view of η(x) being a polynomial), and thus we can always find ϕ − and ϕ + such that ϕ − ≤ ϕ(x) ≤ ϕ + on [0,1], and similarly for Ã(x) ≤ Ã+ .The first term on the RHS of ( 18) can be viewed as a simple e Ã(X T )− Ã+ -coin flip, whilst the second term is precisely the probability that all points in a unit rate Poisson point process Φ we can thin Φ to a Poisson point process on [0, T ] × [0, ϕ + − ϕ − ] and hence simulate an event with probability given by the RHS of ( 18).
This allows for exact paths from the non-neutral Wright-Fisher process to be returned by first simulating the appropriate Poisson point process, subsequently generating draws from the neutral Wright-Fisher process at the time-stamps returned by the Poisson point process, checking whether the generated points all lie in the appropriate region, and and finally running a simple e Ã(X T )− Ã+ -coin flip.
In order to calculate Ã+ , ϕ − and ϕ + , a Polynomial class (with associated root finding algorithm implementing the Jenkins-Traub algorithm, developed by Bill Hallahan1 ) was used.Whilst the implementation of this routine should work for polynomials of any degree, only polynomials η(x) of degree at most 25 were allowed to ensure that the code returns reliable output within a reasonable amount of time.

Monotonic upper and lower bounds for the new denominators
In this section we show that the denominators in the transition densities for both the diffusion (equations ( 8) and ( 10)) and the diffusion bridge (equations ( 12), ( 16) and ( 17), as well as equations (PA1) through to (PC3)) allow for monotonic sequences of upper and lower bounds.By comparing ( 8) and ( 10), as well as ( 12), ( 16) and ( 17), it becomes clear that we can consider solely the denominator ) as the proofs for the other quantities follow using near identical arguments.
We further emphasise once more (as done in Section 3), that for the bridge case we need only need consider the cases (PA1), (PA2), (PA3), and (PC2) in order to be able to simulate draws from any Wright-Fisher diffusion bridge process.Additionally, observe that the denominator for (PA3) is given by q θ 0 (t) for which monotonic upper and lower bounds are immediate, whereas (PC2) is precisely the case covered by Proposition 3 in Jenkins and Spanò (2017).It therefore remains to find monotonically converging sequences of upper and lower bounds for each of: Further, by equation ( 5) in Griffiths et al. (2018), (20) admits the required monotonic bounds through analytic expressions for the falling factorial moments of the ancestral process (see Theorem 5 in Griffiths et al. (2018) and the preceding paragraphs for full details).

Calculations for (21)
Dealing with (21) requires some more work; we start by observing that (1 We can modify the arguments in Lemma 1 in Jenkins and Spanò (2017) to deduce that for L m ∼ Bin(m, x) we have that To see this, observe that for l ≤ ⌊mx⌋ where in the last inequality we used the fact that l ≤ ⌊mx⌋ ≤ mx.When l ≥ ⌊mx⌋, we have that by observing that when mx > 1, m+1 l+1 ≤ m+1 mx ≤ 1+ 1 x , whereas for mx ≤ 1, m+1 l+1 ≤ m+1 ≤ 1+ 1 x .Summing together ( 25) and ( 26

Calculations for (23)
Note first that and observe that the terms inside the inner summation (excluding the factor (−1) k−m ) correspond to the terms b (t,θ) k (m) as defined in Proposition 1 in Jenkins and Spanò (2017).Let , and observe that we can re-write the above as ∞ i=0 (−1) i d i with For ε ∈ (0, 1) fixed, denote by Proposition 3 in Jenkins and Spanò ( 2017) can be restated for the case we consider here as follows.
Proof.The proof proceeds as in Jenkins and Spanò (2017).As m > E t , 2j ≥ C t m−j , and thus by Proposition 1 in Jenkins and Spanò (2017) and summing over j gives The above reasoning also leads to The result follows.
6 Approximations for small times and implementation Whenever the simulation time increments become too small, numerical instabilities crop up when computing contributions to the quantities q θ m (t).Thus (as done in Jenkins and Spanò (2017)), adequate approximations are necessary which make use of the small time asymptotics of q θ m (t).Theorem 4 in Griffiths (1984) gives that as t → 0, the ancestral block counting process of the coalescent is well approximated by a Gaussian random variable with mean and variance (note that Theorem 4 in Griffiths (1984) is missing a factor of β −2 ).In light of this, whenever the time increment t falls below a specific threshold ε G , EWF makes use of the above Gaussian approximation, such that the probabilities q θ m (t) are replaced by their (suitably rounded) Gaussian counterparts.In the current implementation of EWF, the threshold ε G was set to 0.08 after extensive testing as it was found that such a cutoff ensured a suitable trade-off between retaining precision by employing the approximation only when necessary, and having a robust and efficient implementation.
For the diffusion bridge case we apply similar approximations to both q θ m (s) and q θ k (t − s), but we also introduce an additional threshold ε D < ε G below which we approximate draws for the law of a diffusion bridge through draws from the law of a diffusion.This is necessary due to the fact that the mean µ given above for the Gaussian approximation is inversely proportional to the time increment t.Thus if either of the time increments s or t − s is small, the pmf {p m,k,l,j } m,k,l,j∈N spreads out very thinly over N 4 leading to a loss of precision due to the small quantities involved coupled with infeasible run times, even when the above illustrated Gaussian approximations are used.
In such cases (i.e.s < ε D or t − s < ε D ), EWF first simulates a draw from the corresponding Wright-Fisher diffusion started at x and sampled at time s, computes the increment between the generated draw Y ′ and the start point x, and superimposes it onto a linear interpolation between the left and right end-points x and z to generate the required draw Y .The linear interpolation employed explicitly make use of time increments s and t − s to account for the fact that the returned draw Y should come from a diffusion bridge starting at x and ending at z, with appropriate mechanisms in place to ensure that the output remains within the interval , the above detailed (rounded) Gaussian approximations are used for the corresponding {q θ i } i∈N within the appropriate time interval, whilst the standard sampling scheme is used for time increments which exceed ε G .A threshold of 0.008 was chosen for ε D following extensive testing, such that the resulting implementation of EWF retained robustness and efficiency and refrained from using such approximations unless their absence led to unfeasible run times.We mention that both thresholds can be altered if desired through the fields g1984 (for ε G ) and bridgethreshold (for ε D ) of the Options class found in the myHelpers.hheader file (although we would advise against this).

Output validation
Output was validated by generating 10,000 samples for a wide variety of cases and subsequently comparing this to a truncation of the transition density by means of Kolmogorov-Smirnov test as well as QQ-plots.We point out that we present only neutral output here as the non-neutral output is generated using the same rejection procedure as used in Jenkins and Spanò (2017).
To illustrate how the transition density was truncated, consider the case (PC2), which involves a sum over four indices, two of which are infinite.By using an iterative scheme, the mode over these four indices was found and its contribution to the density for a given point y ∈ [0, 1] was calculated.Subsequently the denominator of (PC2) was evaluated up to machine precision, and an appropriate truncation level was chosen by multiplying together the resulting denominator, the mode's contribution to the density and a tolerance parameter.Similar truncations were employed for all the other diffusion and diffusion bridge cases.

Diffusions conditioned on non-absorption
Samples were generated using 9 different parameters setups featuring starting points x ∈ {0, 0.5, 1}, sampling times t ∈ {0.01, 0.05, 0.5}, and mutation parameter θ = (0, 1).The output is plotted below, starting with the case when x = 0, with the sampling time increment t increasing when going left to right across plots.All of the Kolmogorov-Smirnov tests and QQ-plots below confirm that the output is indeed coming from the correct distribution.

Unconditioned diffusions
In the case when the diffusion was allowed to be absorbed at the boundaries, simulations for the following cases were obtained: start points x ∈ {0.25, 0.5, 0.75}, sampling times t ∈ {0.05, 0.25, 0.5}, and mutation parameter θ = 0. We report the probability of being absorbed at either boundary in the table below, where P denotes the empirical estimate for this quantity whereas P is the theoretical value obtained by evaluating the truncation to the transition density at the boundary.All of the estimated probabilities match their theoretical counterparts, and further both the QQ-plots and Kolmogorov-Smirnov tests confirm that the generated draws are coming from the correct distribution.Figure 5: (Top row): Histograms for 10,000 samples generated from the law of a Wright-Fisher diffusion started at x = 0.25 at time 0, sampled at times t = 0.05, 0.25, 0.5 respectively, with the process allowed to be absorbed at the boundaries.Note that samples equal to 0 or 1 are not included in the above histograms, but their relative frequency can be found from the empirical probabilities found in Figure 6: (Top row): Histograms for 10,000 samples generated from the law of a Wright-Fisher diffusion started at x = 0.5 at time 0, sampled at times t = 0.05, 0.25, 0.5 respectively, with the process allowed to be absorbed at the boundaries.Note that samples equal to 0 or 1 are not included in the above histograms, but their relative frequency can be found from the empirical probabilities found in Figure 7: (Top row): Histograms for 10,000 samples generated from the law of a Wright-Fisher diffusion started at x = 0.75 at time 0, sampled at times t = 0.05, 0.25, 0.5 respectively, with the process allowed to be absorbed at the boundaries.Note that samples equal to 0 or 1 are not included in the above histograms, but their relative frequency can be found from the empirical probabilities found in Table 3.The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.
We further considered the following sampling times for each bridge:  4 above, sampled at the times given by the corresponding row in Table 5.The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.

Unconditioned bridges
When the diffusion bridge is allowed to be absorbed at the boundary and θ = 0, we need only consider the cases when z ∈ {0, 1}.To this end we considered the following two setups: Bridge 1 (0,0.25)(0.3,1) Bridge 2 (0,0.5)(0.5,0) Table 6: The left and right endpoints for the three different bridges simulated, where (t 0 , x 0 ) denotes the bridge's start time t 0 and start point x 0 , (t 1 , x 1 ) denotes the second observation time and point for the diffusion bridge and so on.
We further considered the following sampling times: Histograms for 10,000 samples generated from the law of the Wright-Fisher diffusion bridge 'Bridge 1' (allowed to be absorbed at 1) as given in Table 6, sampled at the times given by the corresponding row in Table 7.Note that the samples equal to 1 are not included in the above plots, but their relative frequency can be found in Table 8.The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.6, sampled at the times given by the corresponding row in Table 7.Note that the samples equal to 0 are not included in the above plots, but their relative frequency can be found in Table 9.The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.

Non-neutral diffusions and diffusion bridges
Non-neutral Wright-Fisher paths can be generated (as described in Section 4) through the use of neutral paths coupled with an appropriate Poisson point process.This technique was proposed in Jenkins and Spanò (2017) and is used (without any alteration) in the current implementation of EWF to return non-neutral draws from the laws of both diffusions and diffusion bridges.
Thus, although EWF does allow for non-neutral draws under a very broad class of selective regimes (and instructions on how to do this can be found in the respective configuration files), we omit the resulting output.

Figure 1 :
Figure 1: Illustration of 30 candidate trajectories for the horse coat color data found in Ludwig et al. (2009) simulated using EWF (note that the observed frequencies (black crosses) are assumed to be exact observations of the underlying diffusion).Simulations used the inferred selection coefficient s = 0.0007 with a consensus effective population size N e = 10, 000(Ludwig et al., 2009;Malaspinas et al., 2012;Schraiber et al., 2016), giving σ = 2N e s = 14.We used θ = 0 and a generation time of 5 years.

Figure 2 :Figure 3 :Figure 4 :
Figure2: (Top row): Histograms for 10,000 samples generated from the law of a Wright-Fisher diffusion conditioned on non-absorption, started at x = 0 at time 0, sampled at times t = 0.01, 0.05, 0.5 respectively.The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.

Figure 8 :
Figure 8: (Top row): Histograms for 10,000 samples generated from the law of the Wright-Fisher diffusion bridge 'Bridge 1' in Table4above, sampled at the times given by the corresponding row in Table5.The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.

Figure 11 :
Figure11: (Top row): Histograms for 10,000 samples generated from the law of the Wright-Fisher diffusion bridge 'Bridge 1' (allowed to be absorbed at 1) as given in Table6, sampled at the times given by the corresponding row in Table7.Note that the samples equal to 1 are not included in the above plots, but their relative frequency can be found in Table8.The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.

Figure 12
Figure12: (Top row): Histograms for 10,000 samples generated from the law of the Wright-Fisher diffusion bridge 'Bridge 2' (allowed to be absorbed at 0) as given in Table6, sampled at the times given by the corresponding row in Table7.Note that the samples equal to 0 are not included in the above plots, but their relative frequency can be found in Table9.The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.

Table 1 :
Empirical ( P) and theoretical (P) absorption probabilities for the diffusion started at x = 0.25.

Table 1 .
The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.

Table 2
. The truncated transition density is plotted in red.(Bottom row): QQ-plots for the corresponding samples with the p-value returned from the Kolmogorov-Smirnov test reported above the plot.