Small-time approximation of the transition density for diffusions with singularities. Application to the Wright-Fisher model.
Abstract
The Wright-Fisher (W-F) diffusion model serves as a foundational framework for interpreting population evolution through allele frequency dynamics over time. Despite the known transition probability between consecutive generations, an exact analytical expression for the transition density at arbitrary time intervals remains elusive. Commonly utilized distributions such as Gaussian or Beta inadequately address the fixation issue at extreme allele frequencies (0 or 1), particularly for short periods. In this study, we introduce two alternative parametric functions, namely the Asymptotic Expansion (AE) and the Gaussian approximation (GaussA), derived through probabilistic methodologies, aiming to better approximate this density.
The AE function provides a suitable density for allele frequency distributions, encompassing extreme values within the interval [0,1]. Additionally, we outline the range of validity for the GaussA approximation. While our primary focus is on W-F diffusion, we demonstrate how our findings extend to other diffusion models featuring singularities. Through simulations of allele frequencies under a W-F process and employing a recently developed adaptive density estimation method, we conduct a comparative analysis to assess the fit of the proposed densities against the Beta and Gaussian distributions.
Key words: diffusions with singularities,
Wright-Fisher model,
transition density.
Classification:
1 Introduction
The Wright-Fisher (W-F) model, initially proposed by Wright1931, offers a robust mathematical framework for studying the dynamics of allele frequencies within a population over successive generations.
To begin, we introduce the discrete W-F model alongside its diffusion approximation. The discrete W-F model, as delineated by Wright1931, operates within a population of constant size and two allelic types . It is represented as a Markov chain with state space , where denotes the count of allele at generation . Transition probabilities within this model are defined by a conditional Bernoulli law:
Scaling the number of generations relative to the population size permits approximate the model with a diffusion process. Let represent the Markov chain, emphasizing the population size, with as . Introducing the continuous-time process yields as tends to infinity, being this convergence in distribution, as demonstrated in the book Ethier1986. The limiting process is a Markov diffusion with state space , governed by the following Stochastic Differential Equation (SDE):
(1.0.1) | |||||
Here, denotes a standard Wiener process.
The infinitesimal generator of this diffusion process is defined as:
where is a function twice continuously differentiable and defined on with appropriate boundary conditions.
Let us introduce the main focus of our work. Since no simple analytical expression is available for the distribution of allele frequencies in a given generation, the application of the discrete W-F model in inference has relied on approximate results. An overview of these results can be found in Tataru2016, which compares several approximations of the W-F distribution used in statistical inference under different assumptions.
The Gaussian distribution is commonly employed ( Bonhomme2010; Coop2010; Gautier2010; Pickrell2012; Lacerda2014; Terhorst2015), albeit with limitations due to the discrete nature of the model and the appearance of fixation or loss events at the extremes of the frequency spectrum. Most of the authors cite Nicholson2002 to justify the use of this distribution, but no proof is given in that article. This approximation is not always appropriate due to (i) the discrete nature of the model, (ii) the range of frequencies, which is restricted to the interval , and (iii) the appearance of atoms at 0 and 1 in the allele frequency, corresponding to fixation and loss of an allele, respectively. As the number of generations increases, given a constant population size, the atoms have an increasing probability of remaining fixed at 0 or 1. One possible solution, explored in the work of Tataru2015, is to use a Beta distribution, Beta with spikes at the extremes of the interval , or a truncated Gaussian distribution. The most promising proposal in terms of goodness of fit of the allele frequency distribution is to consider the diffusion model (1.0.1), but since its transition density does not have a closed analytical expression, it is not possible to perform statistical inference. These drawbacks make finding analytical expressions and/or valid ranges of the Gaussian distribution very important for statistical inference.
In this paper we propose a more general framework that allows obtaining a finite definite density function on the interval . Our methodology, inspired by seminal works in the field, provides an asymptotic expansion of the probability density function solution for the W-F diffusion model, particularly effective when is small. To find such a density, we generalize to a singular diffusion the technique developed in the work of DacunhaCastelle1986 and Voronka1975. Other works that complement the scope of the cited work of Voronka and Keller are those of antonelli1978, tier1978, and tier1981, where the Ohta–Kimura model with two-locus di-allelic systems with linkage are studied, for example.
As well as Voronka1975, our method provides an asymptotic expansion of the solution of the probability density function of the forward Kolmogorov equation for the W-F diffusion model that is valid when is small (where is the number of generations and is the constant population size). It is important to point out that our method is probabilistic in essence, instead of asymptotic methods in partial differential equations as is the case of the aforementioned work.
While our primary focus remains on the W-F diffusion, we present a generalized solution applicable to diffusions with singularities parameterized by two parameters, and . When these parameters are equal to , we obtain W-F diffusion. Thus, we show how to develop a general procedure to find the transition density of a diffusion that is state dependent and singular.
To evaluate the efficacy of our proposed density and compare it with existing methods, we conduct simulations of discrete W-F trajectories. For a fixed time , we use an adaptive estimation method to obtain a continuous empirical density from the simulated data (see, e.g., Bertin2011). Having a continuous density from the simulated data allows us to better evaluate the fit of different continuous densities to the data using the Hellinger distance and the norm.
Our work is organized into various sections. In Section 2, we introduce our notation, the underlying model that we are investigating, the process of transforming state dependent diffusion, significant findings in discovering transition densities of diffusion processes, and the lemma that provides the transition density. We provide a comprehensive proof of this lemma in Section 3. For Section 4, we present significant outcomes associated with the W-F model, as well as corresponding proofs and special scenarios, such as neutral, with mutation and selection. Proofs for these scenarios can be found in Section LABEL:apend. Additionally, in Section LABEL:section:_model_evaluation, we assess the suitability of the defined models and two other models commonly used: the Beta and Gaussian distribution.
2 Preliminaries
In this section we begin by defining a stochastic differential equation (SDE) that has a drift component, , and a local variance, , depending on . We then define the adaptation of Lemma 2, presented by DacunhaCastelle1986 paper, so that the transition density can be studied. In this lemma, both the drift term and the local variance term require certain regularity of these functions. Using the Lamperti transform (see, for example, moller2010), we study a transformed process and the conditions for the drift term of this process to fulfill the conditions of the proposed lemma.
We will find in subsection 3.3 (proposition 6) an exact formula for the transition density of the following diffusion model
(2.0.1) |
where is a standard Wiener process. The W-F diffusion is a particular case of this model, when .
As local variance term of this SDE is not constant, so it is necessary to transform it.
2.1 Local variance transformation for a SDE
Through the Lamperti transformation we can convert a SDE, with drift and local variance term dependent on , into a SDE with a constant local variance term.
Let us consider the following stochastic differential equation,
(2.1.1) |
where its drift component is s smooth function, and the local varainace is twice differentiable in , the continuity of this function implies that it is lower bounded by a constant , and ia a standard Wiener process.
Let be a transformation of , where is a twice differentiable function; Itô lemma yields
Defining as
(2.1.2) | ||||
Then the result for the transformed process is
(2.1.3) |
The process is now the solution of a SDE with a constant local variance term equal to one. With the additional hypothesis that has an inverse, i.e , we can re-write the equation (2.1.3) as
(2.1.4) |
Now, we have a SDE, resulting from the Lamperti transformation, with constant local variance equal to one and drift term equal to
In the particular case when we obtain
(2.1.5) |
2.2 Transition density for a diffusion with constant local variance term
The above transformation allows us to work, without loss of generality, with a diffusion that has a constant term of local variance. Now based on Lemma 2, of DacunhaCastelle1986, we present the following lemma that allows obtaining a exact formula for the transition density.
Lemma 1.
[Transition density] Let us consider the following stochastic differential equation
(2.2.1) |
where is the standard Wiener process and is a smooth drift function and verifies
(2.2.2) |
Let be the transition density of the process defined in (2.2.1), and be a standard Brownian bridge and let us denote . Then, can be written as
(2.2.3) |
where .
A detailed proof of Lemma 1 can be found in the following section.
Furthermore, using the above lemma we can obtain an useful approximation for the transition density when .
Lemma 2.
For and , we have
Proof.
Using the approximation when , we get that
Then
(2.2.4) |
∎
Corollary 3.
[Asymptotic expansion] The transition density, , can be written as follows
(2.2.5) |
Above, we have used the Landau’s symbols and , we say the or , if or a constant, respectively.
The equation defined in (2.2.5) gives us an approximation to the transition density that is valid when approaches 0. For this reason, we will refer to this approximation as the Asymptotic Expansion (AE).
Given that we are going to work with a transformed process , it is necessary to establish a link between the results for the original process, and the transformed one. The following lemma is presented in order to state the relation between their two transition densities.
Lemma 4.
Let be a real-valued Markov process with transition density and an increasing and differentiable function. Let us define the process with transition density . Then, for any real numbers , the transition density of can be written as
(2.2.6) |
Proof.
Let be real numbers such that . Then,
{IEEEeqnarray*}rCl
P^x_0(X(t)∈[c_1, c_2])=P^y_0(Y(t)∈[F(c_1),F(c_2)])&=∫_F(c_1)^F(c_2)p^Y_t(y_0,y)dy
=∫_c_1^c_2 p^Y_t(F(x_0),F(x))F’(x)dx,
where and .
By the definition of transition density, we conclude that (2.2.6) holds.
∎
This lemma implies that if we have an exact or approximate analytical expression for the transition density of , we also have one for .
3 Proof for transition density lemma 1
To study the transition density defined in (2.2.3), we further analyze the terms involved in it. Starting with the associated semigroup and how it is related to the drift function, continuing with the associated drift function and finally and analysis of the function considered.
For ease of notation let us write the transition density of as
(3.0.1) |
where
By equation (2.2.6), we obtain
(3.0.2) |
3.1 Semigroup associated
The semigroup associated to the Markov process defined in (2.2.1) is given by
(3.1.1) |
where is a smooth function with compact support and is the expectation conditioning on . The semigroup in (3.1.1) can be written using the Cameron-Martin-Girsanov formula (see Oksendal2000 pp. 123) (this formula holds under the condition when ) as {IEEEeqnarray}rCl P_tf(x) & = E^x[e^∫_0^t μ(W(s))dW(s)-12∫_0^t μ^2(W(s))dsf(W(t))]
Which, when it is applied to (2.1.5), gives
(3.1.2) |
Now, we analyze the drift function involved.
3.2 Drift function
In our case, the drift function is defined by . Considering that and is an increasing function, we have
Considering these results and by means of Itô lemma,
Implying that,
Replacing this result in (3.1.2), it results a formula without terms involving stochastic integrals
(3.2.1) |
Remark 5.
It is important to point out that all the results obtained above can be obtained without difficulty from the paper of DacunhaCastelle1986.
3.3 Analysis of
This subsection aims to extend the precedent results to more general local variance terms. First let us assume that is concave and let us define
The above result enables us to write,
By duality, it results that for , yields,
In addition, for
This relationship means that
(3.3.1) |
Below we extend the class of functions for which the above result holds, based on the following assumptions
-
(A1)
Let us consider that , is concave, and it tends to zero at and , such that .
-
(A2)
We can approximate and its two derivatives by a sequence of functions with two continuous derivatives such that , , as well as its derivatives
These assumptions allow us to approximate the expression (3.2.1) by a similar expression written for the sequence instead of ; then, by the dominated convergence theorem,we can exchange limit with expectation.
To clarify this point, let be defined as , but replacing instead of , and let us denote the transition density. The Kolmogorov differential equation, implies that .
Additionally, we have
(3.3.2) |
where, by equality in distribution, we have substituted, the sequence of Brownian bridges, , by a fixed Brownian bridge .
Considering that is a positive function, we use dominated convergence theorem to show that the term (3.3.2) converges to
Thus, we obtain
and also the relationship (3.3.1), holds.
These computations allow us to give the following result
Proposition 6.
Based on assumptions (A1) and (A2), we get
Remark 7.
We think this is a new result; we have not found a similar one in the literature.. The class of functions defined by , for satisfies the conditions (A1) and (A2). When it corresponds to the W-F diffusion.
4 Wright - Fisher diffusion: Asymptotic Expansion and Gaussian approximation
In the following we concentrate in the main subject of this paper: the W-F diffusion ( ).
A parallel development can be carried out for other values of the parameters, although we will not consider this problem in the present work.
4.1 Analytical expression for the transition density
The goal of this section will be to prove an exact analytical expression for the transition density and its asymptotic expansion for small , of the W-F diffusion. When this diffusion is considered, it produces the following
(4.1.1) | |||
(4.1.2) |
Let us recall that in this case the following function writes . The main result is
Theorem 8.
The transition density of the W-F diffusion can be written as
(4.1.3) |
and its asymptotical expansion as
(4.1.4) |
for all .
Proof.
Let us start by proving . For ease of notation let us write the transition density of , based on Proposition 6, as
(4.1.5) |
where
By equation (2.2.6), we obtain
(4.1.6) |
Let us further develop the terms in equation (4.1.5). For every pair of positive real numbers and such that , we have that
{IEEEeqnarray*}rCl
M(y)-M(y_0) &= -12∫_y_0^ycot(u) du
= 12 (log—siny—- log—siny_0— )
= 14 [logsin^2(y_0)-logsin^2(y) ]
= 14[log