# A Faster and More Accurate Algorithm for Calculating Population Genetics Statistics Requiring Sums of Stirling Numbers of the First Kind

^{*}Infectious Diseases Programme, Department of Medicine, Yong Loo Lin School of Medicine, National University of Singapore, Singapore 119228, Singapore & Infectious Diseases Group, Genome Institute of Singapore, Singapore 138672^{†}Centrum Wiskunde & Informatica (CWI), Science Park 123, 1098 XG Amsterdam, The Netherlands

- 2Corresponding author: Division of Infectious Diseases, Department of Medicine, Yong Loo Lin School of Medicine, National University of Singapore, Singapore 119228 & Laboratory of Bacterial Genomics, Genome Institute of Singapore, Singapore 138672. E-mail: slchen{at}gis.a-star.edu.sg

## Abstract

Ewen’s sampling formula is a foundational theoretical result that connects probability and number theory with molecular genetics and molecular evolution; it was the analytical result required for testing the neutral theory of evolution, and has since been directly or indirectly utilized in a number of population genetics statistics. Ewen’s sampling formula, in turn, is deeply connected to Stirling numbers of the first kind. 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 . By reducing the calculation from a sum of terms involving Stirling numbers to a single estimate, we simultaneously improve accuracy and dramatically increase speed.

- Population genetics statistics
- Evolutionary inference
- Stirling numbers of the first kind
- Asymptotic analysis
- Numerical algorithms
- Cumulative distribution functions

The dominant paradigm in population genetics is based on a comparison of observed data with parameters derived from a theoretical model (Casillas and Barbadilla 2017; Nielsen 2001). 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 of 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 (Nielsen 2001). Such selective pressures are relevant for understanding and modeling practical problems such as influenza evolution over time (Grenfell *et al.* 2004) and during vaccine production (Chen *et al.* 2019); adaptations in human populations, which may impact disease risk (Wollstein and Stephan 2015; Quintana-Murci 2016); and the emergence of new infectious diseases and outbreaks (Wu *et al.* 2016).

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 (Nielsen 2001; Casillas and Barbadilla 2017; Fu 1997). 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 1972), which describes the expected distribution of the number of alleles in a sample of individuals, given the nucleotide diversity.

Ewens’s sampling formula not only was a seminal result for population genetics, but also established connections with combinatorial stochastic processes, algebra, and number theory (Crane 2016). For population genetics, in particular, Ewens’s sampling formula provided a key analytical result that finally enabled mathematical tests of the neutral theory of evolution (Crane 2016; Nielsen 2001). It has given rise to several classical population genetics tests for neutrality, including the Ewens-Watterson test, Slatkin’s exact test, Strobeck’s *S*, and Fu’s (Fu 1997; Strobeck 1987). 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, Fu’s 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 (Wu *et al.* 2016). However, Stirling numbers rapidly grow large and thus explicit calculation can easily overwhelm the standard floating point range of modern computers.

In previous work, an asymptotic estimator for individual Stirling numbers was used to solve the problem of computing Fu’s for large datasets, which are now becoming common due to rapid progress in DNA sequencing technology (Chen 2019). Without such improved numerical methods, Fu’s calculations for data sets as small as 170 sequences can cause overflow, preventing the use of these statistics for genome-wide screens of selection. This algorithm based on estimating individual Stirling numbers solved problems of numerical overflow and underflow, maintained good accuracy, and substantially increased speed compared with other existing software packages (Chen 2019). However, there was still a need to estimate multiple Stirling numbers (up to half the total number of sequences). Here, we explore the potential for further increasing both accuracy and speed in calculating Fu’s by using a single estimator for the entire sum, which involves multiple Stirling numbers.

## Methods

### General definitions and theory

We take a population of *n* individuals, each of which carries a particular DNA sequence (referred to as the allele of individual *i*). 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:(1)We also define a set of unique alleles which have the property of . The ordinality of is denoted *m*, *i.e.*, the number of distinct alleles in the data set.

Building upon on Ewens’s sampling formula (Fu 1997; Ewens 1972), it has been shown that the probability that, for given *n* and *θ*, at least *m* alleles would be found, is(2)where is the Pochhammer symbol, defined by(3) is a Stirling number and is defined by:(4)Fu’s is then defined as:(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 (Chen 2019).

Because of the relation in (4), the statistics quantity satisfies . Also, this relation and (3) show that are non-negative. We have the special values(6)There is a recurrence relation(7)which easily follows from (4). For a concise overview of properties, with a summary of the uniform approximations, see (Gil *et al.* 2007, §11.3).

We introduce a complementary relation(8)leading to an alternate calculation for Fu’s of(9)The recent algorithm considered in (Chen 2019) is based on asymptotic estimates of derived in (Temme 1993), which are valid for large values of *n*, with unrestricted values of . It avoids the use of the recursion relation given in (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 *n* and *m*: to find *θ* either from the equation , when is given, or from the equation , when is given.

### Remarks on computing

When computing the quantity defined in (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 6.6561 when using the first relation in (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 (8) and using the second relation in (9). A simple criterion to decide about this can be based on using the saddle point (see Remark 1 below).

A second point is numerical overflow when *n* is large, because rapidly becomes large when *m* is small with respect to *n*. For example, when , we have(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) will also be large with *n*; we have .

We can write the sum in (2) in the form(11)Leading to a corresponding modification in the recurrence relation in (7) for the scaled Stirling numbers:(12)To control overflow, we can consider the ratio(13)This function satisfies if . For small values of *n* we can use recursion in the form(14)For large values of *n* and all , we can use a representation based on asymptotic forms of the gamma function.

It should be observed that using the recursions in (7) and (12) is a rather tedious process when *n* is large. For example, when we use it to obtain for all , we need all previous with for all . A table look-up for in floating point form may be a solution. When *n* is large enough, the algorithm mentioned in (Chen 2019) evaluates each needed Stirling number by using the asymptotic approximation derived in (Temme 1993).

### Data availability

Code implementing the new estimator for Fu’s in R is available at https://github.com/swainechen/hfufs.

## Results And Discussion

### Analytical results

The new algorithm is based on the following results, which we describe in two theorems.

**Theorem 1.** The statistics quantity introduced in (2) has the representation as an integral in the complex z-plane(15)where *n* and *m* are positive integers, , *θ* is a real positive number, and is a circle at the origin with radius . The symbol denotes the Pochhammer symbol introduced in (3).

Observe that we have raised in the parameters *n* and *m* with unity; this is convenient in the mathematical analysis. The proof of this theorem will be given in the Appendix (Proof of Theorem 1).

**Corollary 1.** The complementary quantity introduced in (8) has the representation

The main asymptotic result is given in the second theorem.

**Theorem 2.** *has the representation*(17)where is the incomplete beta function defined by(18)with(19)The term is a function of which we give a one-term approximation in (32).

**Corollary 2.** The complementary quantity has the representation

This follows from Theorem 1 and the complementary relation of the incomplete beta function

(21)Note also that the incomplete beta function in (17) has the representation (see (Paris, 2010, §8.17(i)))(22)and from the complementary relation in (21) it follows that the function in (20) has the expansion

(23)The representation in this theorem in terms of the probability function shows the characteristic role of as a cumulative distribution function of the Stirling numbers. The representation can also be viewed as an asymptotic representation in which the incomplete beta function is the main approximant.

The proof of Theorem 2 can be found in the Appendix (Proof of Theorem 2), but we give here some preliminary information about functions used in the proof to explain the definition of the parameter *τ* in (17). It is a function of *θ* and arises in certain transformations of the integral given in Theorem 1. For this we need the function(24)and its derivative(25)With the function we can write (15) in the form(26)Then the saddle point of the integral in (26) follows from the equation(27)There is a positive saddle point when .

Next to these functions we introduce a function for complex values of a variable *t*:(28)where . These functions are related by(29)with condition . In this way, using this relation as a transformation of the variable *z* to *t*, we can write (26) as(30)The parameter *τ* in Theorem 2 is defined as the positive solution of the equation(31)In Figure 1 we show the graphs of (Figure 1A) and (Figure 1B) for , . For these values the saddle points are and . The sign condition for the relation in (29) means that the left branches of the convex curves correspond with functions values for and , and the right branches with values for and . Clearly, we have a one-to-one relation between the positive *z* and *t*-variables.

A first-order approximation of the function in (17) and (20) reads(32)where(33)and the function is defined in (30). The value follows from evaluating (see (30)) at , by observing that both functions and vanish when (hence, ). Then, l’Hôpital’s rule can be used to obtain .

In Figure 2 we show the error curves in (34) for Fu’s (9) for . We show examples for , (Figure 2A) and , (Figure 2B). The solid curves are for when using , the dashed curves when using with the asymptotic estimate given in (32). For ease of visualization, the error has been multiplied by a factor 10 or 100 in Figure 2. We have used the following *mollified error* function(34)where is the approximation of . This mollified error is exactly the relative error unless is small. Because will vanish when (which also means that *θ* is near the transition value (in Figure 2A) and (in Figure 2B) (see Remark 1)), we cannot use relative error for all . This explains the non-smooth curves in Figure 2.

The final estimator is based on the representations in (17) and (20) and the first order approximation in (32), which are used to calculate Fu’s with one of the two relations in (9) depending on whether , decided as described above.

### Implementation and numerical results

We first summarize the steps to compute Fu’s by using (9) and the first-order approximations (see (32) and (17) or (20))(35)or(36)for large *n*, and .

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

2 With , the positive zero of (see (28)), compute

*τ*, the solution of the equation (see (31))(37)with defined in (24) and defined in (28). 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 (35), and from the first relation in (9).

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

In Table 1, we show the relative errors in the computation of defined in (5). The values of *n*, *m*, and *θ* correspond with those in Table 1 of (Chen 2019). The asymptotic result is from (35). Computations were done with Maple, with Digits = 16. The “exact” values were 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 (Chen 2019). 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 (35) was well controlled at for nearly 99% of the calculations; for 411 calculations where the previous hybrid estimator had an error , the estimate in (35) was more accurate in all but one case (; 3.08e-3 relative accuracy using (Chen 2019); 3.32e-3 relative accuracy using (35)) (Figure 3). Further analysis of the relative error demonstrated that it peaks at intermediate values of , depending on *θ*. These correspond to parameter choices near the transition values , where *t* approaches and *z* approaches in the calculation; notably, they remain well controlled (all values mollified error) regardless of *θ*. The asymptotic behavior (lower relative error) can also be seen as both *n* and *m* increase in Figure 3B.

The fewer calculations led to a clear improvement in calculation speed (median 54.6x faster; Figure 4). The speedup also depends on the parameter choices; in general, the speed advantage is greater when the hybrid calculator requires many calculations (namely, when *m* is small relative to *n*, as the hybrid calculator performs the sum in (2)) (Figure 4B).

## 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 *S* have been difficult to calculate on modern, large data sets using previous methods; we now further improve both accuracy and speed for the calculation of Fu’s for such data sets, using the main estimator in (35). Our plan for a paper about the ability to invert the calculation provides additional future directions in understanding the performance of these statistics. Therefore, the methods used herein may be useful for the development of new statistics that more effectively capture different types of selection.

## Acknowledgments

The authors thank the reviewers for their helpful comments and suggestions. SLC acknowledges Shyam Prabhakar and members of the Chen lab for fruitful discussions. NMT acknowledges CWI, Amsterdam, for scientific support; he thanks Edgardo Cheb–Terrab (Maplesoft) for his help with the Maple code. 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).

## APPENDIX

### Proof of Theorem 1

We use the integral representation of the Stirling numbers that follows from the definition given in (4). That is, by using Cauchy’s formula,(38)where is a circle around the origin with radius *R*. We can take *R* as large as we like. As in (Temme, 1993, §3), it is convenient to proceed with(39)Using the definition of in (2) we have(40)and using (39) we obtain(41)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(42)This concludes the proof of Theorem 1.

Corollary 1 now follows by using the theory of integrals of analytical functions on complex contours. We have assumed that , but we can take while picking up the residue at . The result is(43)This gives the relation in Corollary 1.

**Remark 1.** When *θ* crosses the value , becomes (almost) . Especially when the parameters *m* and *n* are large, starts with very small values for small *θ*, becomes close to when , and quickly becomes 1 as *θ* increases. We call the transition value for *θ*.

For fixed values of *n* and *θ*, there is also a transition value for *m*; refer to this transition value as . When *n* is large, starts at values near 1 for small *m*, it becomes about when *m* nears , and it becomes quickly small as .

### Proof of Theorem 2

The relation in (29) between the functions (see (24)) and (see (28)) can be used as a transformation of the variable *z* to *t*, as in (Temme, 1993, §3). The result is the integral representation in (30). In Figure 1 we have shown the relationship between *z* and *t*.

The function in (30) has a pole in the *t*-domain; refer to this pole as . This then corresponds with the pole at in the *z*-domain. The relation between *τ* and *θ* follows from the transformation given in (29). In other words, *τ* is defined by the equation(44)where the sign-convention follows from the one used for (29). We can express the existence of the pole of the function defined in (30) by writing(45)In asymptotic analysis, the presence of such a pole is of great interest, especial when (in the *t*-domain) the saddle point (here ) is close to a pole (here *τ*), or even when these points coalesce. See, for example, (Temme, 2015, Chapter 21). Usually, the error function is introduced to handle the asymptotic analysis; in the present case, we use an incomplete beta function. We split off the pole from and write(46)where we assume that is well defined at . To find *A* we use the analytical relation in (29) between *t* and *z*, in particular at (or ). Applying l’Hôpital’s rule, we find that as , which gives *1*. Hence, substituting this form of in (30), we find(47)The radius of the circle in the first integral is larger than *τ*. For the second integral, we take a circle C around the origin such that the singularities of are outside the circle.

In Proof of the incomplete beta relation below, we prove that the first integral in (47) can be evaluated in terms of the incomplete beta function as shown in Theorem 2. We can then write (47) as(48)where(49)A first-order approximation of (49) follows from replacing by its value at the saddle point . This gives(50)where(51)This expression of follows from (33). In a future publication we will give details about the complete asymptotic expansion of the term .

### Proof of the incomplete beta relation

We give a proof of the claim that the incomplete beta function in (48) equals the first integral in (47). That is,(52)where is a circle at the origin with radius larger than *τ*. We have, using the definition of in (28),(53)which is the relation in (22). In the second line we have used a finite number of terms of the infinite expansion of because terms with do not give a contribution.

## Footnotes

*Communicating editor: A. Kern*

- Received March 9, 2020.
- Accepted August 23, 2020.

- Copyright © 2020 Chen, Temme

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.