    # A faster and more accurate algorithm for calculating population genetics statistics requiring sums of Stirling numbers of the first kind

Stirling numbers of the first kind are used in the derivation of several population genetics statistics, which in turn are useful for testing evolutionary hypotheses directly from DNA sequences. Here, we explore the cumulative distribution function of these Stirling numbers, which enables a single direct estimate of the sum, using representations in terms of the incomplete beta function. This estimator enables an improved method for calculating an asymptotic estimate for one useful statistic, Fu's F_s. By reducing the calculation from a sum of terms involving Stirling numbers to a single estimate, we simultaneously improve accuracy and dramatically increase speed.

## Authors

12/17/2020

### A new asymptotic representation and inversion method for the Student's t distribution

Some special functions are particularly relevant in applied probability ...
10/10/2017

### On some difficulties in the addition of trapezoidal ordered fuzzy numbers

At the first, we revise the Kosinski definition of the sum of ordered fu...
01/10/2022

### An examination of the spillage distribution

We examine a family of discrete probability distributions that describes...
06/30/2020

### A Computational Criterion for the Irrationality of Some Real Numbers

In this paper, we compute the asymptotic average of the decimals of some...
08/31/2020

### Discrete convolution statistic for hypothesis testing

The question of testing for equality in distribution between two linear ...
10/14/2018

### q-Stirling numbers arising from vincular patterns

The distribution of certain Mahonian statistic (called BAST) introduced ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The dominant paradigm in population genetics is based on a comparison of observed data with parameters derived from a theoretical model [1, 2]. Specifically for DNA sequences, many techniques have been developed to test for extreme relationships between average sequence diversity (number of DNA differences between individuals) and the number alleles (distinct DNA sequences in the population). In particular, such methods are widely used to predict selective pressures, where certain mutations confer increased or decreased survival to the next generation . Such selective pressures are relevant for understanding and modeling practical problems such as influenza evolution over time  and during vaccine production ; adaptations in human populations, which may impact disease risk [5, 6]; and the emergence of new infectious diseases and outbreaks .

Many population genetics tests are therefore formulated as unidimensional test statistics, where the pattern of DNA mutations in a sample of individuals is reduced to a single number

[2, 1, 8]

. Such statistics are heavily informed by combinatorial sampling and probability distribution theories, many of which are built upon the foundational Ewens’s sampling formula

. Ewens’s sampling formula describes the expected distribution of the number of alleles in a sample of individuals, given the nucleotide diversity. Calculation of subsets of this distribution are useful for testing deviations of observed data from a null model; such subsets often require the calculation of Stirling numbers of the first kind (hereafter referred to simply as Stirling numbers). In particular, two population genetics statistics, the Fu’s and Strobeck’s statistics, utilize this approach [8, 10]. The former has recently been shown to be potentially useful for detecting genetic loci under selection during population expansions (such as an infectious outbreak) both in theory and in practice . However, Stirling numbers rapidly grow large and overwhelm the standard floating point range of modern computers.

In previous work, an asymptotic estimate for individual Stirling numbers was used to solve the problem of computing Fu’s for large datasets that are now becoming common due to rapid progress in DNA sequencing technology . This new algorithm solved problems of numerical overflow and underflow, maintained good accuracy, and substantially increased speed compared with other existing software packages. However, the estimation of individual Stirling numbers led to the use of an estimator at least and at most times. Here, we explore the potential for further increasing both accuracy and speed in calculating Fu’s by using a single estimator.

The new estimator for Fu’s has been implemented in R and is available at https://github.com/swainechen/hfufs.

## 2 Background Theory

### 2.1 General definitions

We take a population of individuals, each of which carries a particular DNA sequence (referred to as the allele of individual ). We define a metric, to be the number of positions at which sequence differs from . Then, we denote the average pairwise nucleotide difference as (hereafter referred to simply as ), defined as:

 θ=2n(n+1)n−1∑i=1n∑j=i+1dist(Di,Dj). (2.1)

We also define a set of unique alleles which have the property of . The ordinality of is denoted , i.e. the number of distinct alleles in the data set.

Building upon on Ewens’s sampling formula [8, 9], it has been shown that the probability that, for given and , at least alleles would be found, is

 S′n,m(θ)=1(θ)nn∑k=m(−1)n−kS(k)nθk,θ>0, (2.2)

where is the Pochhammer symbol, defined by

 (θ)0=1,(θ)n=θ(θ+1)⋯(θ+n−1)=Γ(θ+n)Γ(θ). (2.3)

is a Stirling number and is defined by:

 (θ)n=n∑k=0(−1)n−kS(k)nθk, (2.4)

Fu’s is then defined as:

 Fs=lnS′n,m(θ)1−S′n,m(θ). (2.5)

Fu’s thus measures the probability of finding a more extreme (equal or higher) number of alleles than actually observed. It requires computing a sum of terms containing Stirling numbers, which rapidly become large and therefore impractical to calculate explicitly even with modern computers .

Because of the relation in (2.4), the statistics quantity satisfies . Also, this relation and (2.3) show that are non-negative. We have the special values

 S(n)n=1 (n≥0),S(0)n=0 (n≥1),S(1)n=(−1)n−1(n−1)! (n≥1). (2.6)

There is a recurrence relation

 S(k)n+1=S(k−1)n−nS(k)n, (2.7)

which easily follows from (2.4). For a concise overview of properties, with a summary of the uniform approximations, see [12, §11.3].

We introduce a complementary relation

 T′n,m(θ)=1−S′n,m(θ)=1(θ)nm−1∑k=0(−1)n−kS(k)nθk, (2.8)

leading to an alternate calculation for Fu’s of

 Fs=lnS′n,m(θ)1−S′n,m(θ)=ln1−T′n,m(θ)T′n,m(θ). (2.9)

The recent algorithm considered in  is based on asymptotic estimates of derived in , which are valid for large values of , with unrestricted values of . It avoids the use of the recursion relation given in (2.7).

In the present paper we derive an integral representation of and of the complementary function , for which we can use the same asymptotic approach as for the Stirling numbers without calculating the Stirling numbers themselves. From the integral representation we also obtain a representation in which the incomplete beta function occurs as the main approximant. In this way we have a convenient representation, which is available as well for many classical cumulative distribution functions. We show numerical tests based on a first-order asymptotic approximation, which includes the incomplete beta function. In a future paper we give more details on the complete asymptotic expansion of , and, in addition, we will consider an inversion problem for large and : to find either from the equation , when is given, or from the equation , when is given.

### 2.2 Remarks on computing S′n,m(θ)

When computing the quantity defined in (2.5), numerical instability may happen when is close to 1. In that case, the computation of suffers from cancellation of digits. For example, take , , . Then , and becomes about when using the first relation in (2.9). However, when we calculate and use the second relation, then we obtain the more reliable result .

We conclude that, when , it is better to switch and obtain from the sum in (2.8), and by using the second relation of in (2.9). A simple criterion to decide about this can be based on using the saddle point (see Remark 3.1 below).

A second point is the overflow in numerical computations when is large, because of the large values of when is small with respect to . For example, when , we have

 S(5)10=−n!(m+5)(m+4)(3m2+23m+38)11520(m−1)!=−269325. (2.10)

Therefore, it is convenient to scale the Stirling number in the form . In addition, the Pochhammer term in front of the sum in (2.2) will also be large with ; we have .

We can write the sum in (2.2) in the form

 S′n,m(θ)=n!(θ)nn∑k=m(−1)n−kˆS(k)nθk,ˆS(k)n=S(k)nn!. (2.11)

Leading to a corresponding modification in the recurrence relation in (2.7) for the scaled Stirling numbers:

 ˆS(m)n+1=1n+1(ˆS(m−1)n−nˆS(m)n). (2.12)

To control overflow, we can consider the ratio

 fn(θ)=n!(θ)n=Γ(n+1)Γ(θ)Γ(θ+n). (2.13)

This function satisfies if . For small values of we can use recursion in the form

 fn+1(θ)=n+1n+θfn(θ),n=0,1,2,…,f0(θ)=1. (2.14)

For large values of and all we can use a representation based on asymptotic forms of the gamma function.

###### Remark 2.1.

It should be observed that using the recursion in (2.7) and (2.12) is a rather tedious process when is large. For example, when we use it to obtain for all , we need all previous with for all . Table look-up for in floating point form may be a solution. When is large enough, the algorithm mentioned in  evaluates each needed Stirling number by using the asymptotic approximation derived in .

## 3 Results and Discussion

### 3.1 An integral representation of S′n,m(θ)

We use the integral representation of the Stirling numbers that follows from the definition given in (2.4). That is, by using Cauchy’s formula,

 (−1)n−mS(m)n=12πi∫CR(z)ndzzm+1, (3.1)

where is a circle around the origin with radius . We can take as large as we like. As in [13, §3], it is convenient to proceed with

 (−1)n−mS(m+1)n+1=12πi∫CR(z+1)ndzzm+1. (3.2)

We derive an integral representation of

 Sn+1,m+1(θ)=1(θ)n+1n+1∑k=m+1(−1)n+1−kS(k)n+1θk=1(θ+1)nn∑k=m(−1)n−kS(k+1)n+1θk. (3.3)

We use (3.2) and obtain

 S′n+1,m+1(θ)=1(θ+1)nn∑k=mθk2πi∫CR(z+1)nzk+1dz. (3.4)

We can take to have on the circle , and we can perform the summation to , because all terms with do not give contributions. In this way we obtain the requested integral representation

 S′n+1,m+1(θ)=θm(θ+1)n12πi∫CR(z+1)nzmdzz−θ,R>θ. (3.5)

To obtain this result we need , but in the integral representation we can take when we pick up the residue at . The result is

 S′n+1,m+1(θ)=1−θm(θ+1)n12πi∫CR(z+1)nzmdzθ−z,R<θ, (3.6)

and we find for (see (2.8))

 (3.7)

For the asymptotic analysis we write (3.5) in the form

 S′n+1,m+1(θ)=e−ϕ(θ)2πi∫CReϕ(z)dzz−θ,R>θ, (3.8)

where

 ϕ(z)=ln((z+1)n)−mlnz=lnΓ(z+n+1)−lnΓ(z+1)−mlnz. (3.9)

Then the saddle point of the integral in (3.8) follows from the equation

 ϕ′(z)=ψ(z+n+1)−ψ(z+1)−mz=0,ψ(z)=Γ′(z)Γ(z). (3.10)

There is a positive saddle point when .

###### Remark 3.1.

When crosses the value , becomes (almost) . Especially when the parameters and are large, starts with very small values for small , its values is about when and it becomes quickly 1 as increases. We call the transition value for .

For fixed values of there is also a transition value for , say, . When is large, starts at values near 1 for small , it becomes about when crosses the transition value , and it becomes quickly small as .

### 3.2 An asymptotic representation of S′n,m(θ)

We use the transformation, as in [13, §3],

 ϕ(z)=χ(t)+B,χ(t)=nln(1+t)−mlnt, (3.11)

with condition , where is the saddle point in the -domain and also the zero of

 χ′(t)=(n−m)t−mt(1+t)=(n−m)t−t0t(1+t). (3.12)

Also

 B=ϕ(z0)−χ(t0),t0=mn−m. (3.13)

With this choice of , the variables and correspond with each other at the respective saddle points.

Using the transformation we obtain

 S′n+1,m+1(θ)=eC2πi∫CS(t+1)ntmf(t)dt, (3.14)

where

 C=B−ϕ(θ)=ϕ(z0)−χ(t0)−ϕ(θ), (3.15)

and

 f(t)=1z−θdzdt,dzdt=χ′(t)ϕ′(z). (3.16)

The contour runs around the origin and includes a pole at that corresponds with the pole in the -plane at .

### 3.3 A representation in terms of the incomplete beta function

The integrands of the integral representations of have a pole at . For the contour integrals this is not a complication, because by using the theory of analytic functions we can deform the contour to avoid the pole, and we can even cross the pole and pick up the residue as we did to obtain the representation in (3.6).

For the integral in the -domain given in (3.14) the same can be done. The function has a pole at a point , say, that follows from the transformation given in (3.11). That means, is defined by the equation

 ϕ(θ)=χ(τ)+B,orϕ(θ)−ϕ(z0)=χ(τ)−χ(t0), (3.17)

and we can show the existence of the pole of the function defined in (3.16) writing

 f(t)=1z−θdzdt=t−τz−θdzdt1t−τ. (3.18)

In asymptotic analysis the presence of such a pole is of great interest, especial when (in the -domain) the saddle point (here ) is close to a pole (here ), or even when these points coalesce. See, for example, [14, Chapter 21]. Usually, the error function is introduced to handle the asymptotic analysis, in the present can we use an incomplete beta function.

We split off the pole from and write

 f(t)=At−τ+g(t), (3.19)

where we assume that is well defined at . To find we use the analytical relation in (3.11) between and , in particular at (or ). Applying l’Hôpital’s rule, we conclude that as , which gives . Hence, substituting this form of in (3.14), we find

 S′n+1,m+1(θ)=e−χ(τ)2πi∫CS(t+1)ntmdtt−τ+e−χ(τ)2πi∫C(t+1)ntmg(t)dt, (3.20)

where we have used (see (3.15) and (3.17))

 C=−ϕ(θ)+ϕ(z0)−χ(t0)=−χ(τ). (3.21)

The radius of the circle in the first integral is larger than , for the second integral we take a circle around the origin such that the singularities of are outside the circle.

The first integral can be evaluated in terms of the incomplete beta function defined by

 Iy(p,q)=1B(p,q)∫y0tp−1(1−t)q−1dt,0

where is the complete beta function

 B(p,q)=Γ(p)Γ(q)Γ(p+q), (3.23)

We will show in the Appendix that

 Iτ1+τ(m,n−m+1)=e−χ(τ)2πi∫CS(t+1)ntmdtt−τ. (3.24)

Hence,

 S′n+1,m+1(θ)=Iτ1+τ(m,n−m+1)+R′n+1,m+1(θ), (3.25)

where

 R′n+1,m+1(θ)=e−χ(τ)2πi∫CS(t+1)ntmg(t)dt. (3.26)

A first-order approximation of this function follows from replacing by its value at the saddle point . This gives

 R′n+1,m+1(θ)∼e−χ(τ)(nm−1) g(t0), (3.27)

where

 g(t0)=f(t0)−1t0−τ,f(t0)=1z0−θ√χ′′(t0)ϕ′′(z0). (3.28)

This expression of follows from the definition of given in (3.16). In a future publication we will give details about the complete asymptotic expansion of the term .

For the complementary function (see (2.8)) we obtain

 T′n+1,m+1(θ)=I11+τ(n−m+1,m)−R′n+1,m+1(θ), (3.29)

where we have used

 Iy(p,q)=1−I1−y(q,p). (3.30)
###### Remark 3.2.

The incomplete beta function in (3.25) has the representation (see [15, §8.17(i)])

 Iτ1+τ(m,n−m+1)=(1+τ)−nn∑j=m(nj)τj, (3.31)

and from the complementary relation in (3.30) it follows that the function in (3.29) has the expansion

 I11+τ(n−m+1,m)=(1+τ)−nm−1∑j=0(nj)τj. (3.32)

### 3.4 Numerical tests

We summarize the steps to obtain the first-order approximations (see (3.25) or (3.29) and (3.27))

 S′n+1,m+1(θ)∼Iτ1+τ(m,n−m+1)+e−χ(τ)(nm−1) g(t0) (3.33)

or

 T′n+1,m+1(θ)∼I11+τ(n−m+1,m)−e−χ(τ)(nm−1) g(t0), (3.34)

for given , and , and to compute Fu’s by using (2.9).

1. Compute the saddle point , the positive zero of ; see (3.10).

2. With , the positive zero of (see (3.12)), compute , the solution of the equation (see (3.17))

 χ(τ)=χ(t0)+ϕ(θ)−ϕ(z0), (3.35)

with defined in (3.9) and defined in (3.11). When there is one solution . When there are two positive solutions, and we take the one that satisfies the condition .

3. When , hence , compute the approximation of by using (3.33), and from the first relation in (2.9).

4. When , hence, , compute the approximation of by using (3.34), and from the second relation in (2.9).

In Table 1 we give the relative errors in the computation of defined in (2.5). The values of , , and correspond with those in Table 1 of . We have used the asymptotic result (3.27). Computations are done with Maple, with Digits = 16. The ”exact” values are obtained by using Maple’s code for , which computes the Stirling numbers of the first kind.

We additionally performed a comparison with the recently published algorithm in . We performed 10,000 calculations with each algorithm and compared the results with an exact calculator. As expected, since the previous algorithm required estimating a Stirling number for each term of the sum, while the current asymptotic estimate directly calculates the sum, both error and compute speed were improved. Relative error for the single term estimate in (3.25) was well controlled at for nearly 99% of the calculations; for 411 calculations where the previous hybrid estimator had an error , the estimate in (3.25) was more accurate in all but one case (; 3.08e-3 relative accuracy using ; 3.32e-3 relative accuracy using (3.25)) (Figure 1). The fewer calculations led to a clear improvement in calculation speed (median 54.6x faster; Figure 2). Figure 1: Comparison of relative error of the estimator from  and the single term asymptotic estimator in (3.25). Relative error for each is calculated against the arbitrary precision implementation described in . In total, 10,000 calculations were performed with nrandomly sampled from a uniform distribution between 50 and 500; m between 2 and n; and θ between 1 and 50. A solid diagonal line is drawn at y=x. Dotted lines are drawn at a relative error of 0.001. Numbers within each quadrant defined by the dotted lines indicate the number of points in each quadrant. The red dot indicates the one case where the relative error was >0.001 and the error of (3.25) was greater than the estimator from . Figure 2: Comparison of run times between the hybrid algorithm from  and the single term asymptotic estimator in (3.25). 100 iterations were run, each with 10,000 calculations. The same set of parameters were used for each algorithm. The order of running the algorithms was alternated with each iteration. The dark horizontal line indicates the median, the box indicates the first and third quartiles, the whiskers are drawn at 1.5x the interquartile range, and outliers are represented by open circles. The median for the hybrid algorithm is 62.64 s; the median for the asymptotic algorithm is 1.17 s.

## 4 Conclusion

The rapid growth of sequencing data has been an enormous boon to population genetics and the study of evolution. Traditional population genetics statistics are still in common use today. The statistics Fu’s and Strobeck’s have been difficult to calculate using previous methods; we now further improve both accuracy and speed for the calculation of Fu’s for large, modern data sets, using the main estimator in (3.25). Our plan for a paper about the ability to invert the calculation provides additional future directions in understanding the performance of these statistics, and the methods used herein may be useful for the development of new statistics that more effectively capture different types of selection.

## Acknowledgments

SLC acknowledges Shyam Prabhakar and members of the Chen lab for fruitful discussions. NMT acknowledges CWI, Amsterdam, for scientific support.
SLC was supported by the National Medical Research Council, Ministry of Health, Singapore (grant numbers NMRC/OFIRG/0009/2016 and
NMRC/CIRG/1467/2017).
NMT was supported by the Ministerio de Ciencia e Innovación, Spain, projects MTM2015-67142-P (MINECO/FEDER, UE) and
PGC2018-098279-B-I00 (MCIU/AEI/FEDER, UE). The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, figures, and tables.

## 5 Appendix

We give a proof of the incomplete beta integral in (3.17). We use the integral representation of the hypergeometric function (see [16, §15.6])

where the contour starts at the origin, encircles the point in the anti-clockwise direction, and returns to the origin. The main conditions are and that is outside the contour.

We also use the relation between the incomplete beta function and the -function (see [15, §8.17(ii)])

 Ix(p,q)=xp(1−x)qpB(p,q)2F1(1,p+q;p+1;x). (5.2)

It follows that

 Ix(p,q)=xp(1−x)q2πi∫(1+)0sp+q−1(s−1)−q(1−xs)ds, (5.3)

and after the substitution and writing , we obtain with , , and as in (3.11)

 I11+τ(p,q)=e−χ(τ)2πi∫C(t+1)ntmdtτ−t, (5.4)

where the pole at is outside the contour. We modify the contour and pick up the residue. In this way we find the relation in (3.17).

## References

•  S. Casillas and A. Barbadilla. Molecular Population Genetics. Genetics, 205(3):1003–1035, Mar 2017.
•  R. Nielsen. Statistical tests of selective neutrality in the age of genomics. Heredity (Edinb), 86(Pt 6):641–647, Jun 2001.
•  B. T. Grenfell, O. G. Pybus, J. R. Gog, J. L. Wood, J. M. Daly, J. A. Mumford, and E. C. Holmes. Unifying the epidemiological and evolutionary dynamics of pathogens. Science, 303(5656):327–332, Jan 2004.
•  H. Chen, J. J. S. Alvarez, S. H. Ng, R. Nielsen, and W. Zhai. Passage Adaptation Correlates With the Reduced Efficacy of the Influenza Vaccine. Clin. Infect. Dis., 69(7):1198–1204, Sep 2019.
•  A. Wollstein and W. Stephan. Inferring positive selection in humans from genomic data. Investig Genet, 6:5, 2015.
•  L. Quintana-Murci. Understanding rare and common diseases in the context of human evolution. Genome Biol., 17(1):225, 11 2016.
•  Z. Wu, B. Periaswamy, O. Sahin, M. Yaeger, P. Plummer, W. Zhai, Z. Shen, L. Dai, S. L. Chen, and Q. Zhang. Point mutations in the major outer membrane protein drive hypervirulence of a rapidly expanding clone of Campylobacter jejuni. Proc. Natl. Acad. Sci. U.S.A., 113(38):10690–10695, 09 2016.
•  Y. X. Fu. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics, 147(2):915–925, Oct 1997.
•  W. J. Ewens. The sampling theory of selectively neutral alleles. Theor Popul Biol, 3(1):87–112, Mar 1972.
•  C. Strobeck. Average number of nucleotide differences in a sample from a single subpopulation: a test for population subdivision. Genetics, 117(1):149–153, Sep 1987.
•  S. L. Chen. Implementation of a Stirling number estimator enables direct calculation of population genetics tests for large sequence datasets. Bioinformatics, 35(15):2668–2670, 2019.
•  A. Gil, J. Segura, and N. M. Temme. Numerical Methods for Special Functions. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2007.
•  N. M. Temme. Asymptotic estimates of Stirling numbers. Stud. Appl. Math., 89(3):233–243, 1993.
•  N. M. Temme. Asymptotic methods for integrals, volume 6 of Series in Analysis. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2015.
•  R. B. Paris. Chapter 8, Incomplete gamma and related functions. In NIST Handbook of Mathematical Functions, pages 173–192. U.S. Dept. Commerce, Washington, DC, 2010.
•  A. B. Olde Daalhuis. Chapter 15, Hypergeometric function. In NIST Handbook of Mathematical Functions, pages 383–401. Cambridge University Press, Cambridge, 2010.