{'fraction_non_alphanumeric': 0.029373092857674466, 'fraction_numerical': 0.012453794437966707, 'mean_word_length': 4.733997155049787, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
What's the value in more moments?

The moment-based non-parametric species richness estimator of Chao is one of the most widely used estimators for the number of unobserved species in a sampling experiment. This is due in large part to its simplicity and robustness. This simplicity can also be a drawback, as it only uses a small amount of information contained in the observed experiment, essentially only the first moment. Previous authors, specifically Harris and Chao, have presented a general moment-based framework for estimating species richness that includes the Chao estimator. The application of this framework has been stymied by both the lack of deep sampling experiments, where higher moments can be accurately estimated, and the lack of efficient algorithms to properly use this information. Technological advances have filled the former void, allowing for sampling experiments orders of magnitude larger than previously considered. We aim to address the latter by connecting results from the theory of moment spaces and Gaussian quadrature to provide a general moment-based non-parametric estimator of species richness that uses more information through more moments and is computationally efficient. We show this estimator performs well and improves upon the Chao estimator on discrete abundance distributions, the simplest cases of heterogeneity. We demonstrate the performance on a simulated populations taken from emerging high-throughput technologies such as RNA-seq, immune repertoire, and metagenomic sequencing.


Introduction
Consider a population composed of an unknown number of distinct species or classes. A sample of individuals is taken with replacement from the population. Once sampled, the species of an individual is known. Following the sampling experiment, the number of sampled individuals from each observed species can be tallied. The problem of estimating the total number of species in the population using information contained in the observed sample is commonly known as the species richness problem. This terminology originated in the context of ecological experiments, but the underlying statistical question is simply to estimate the size of a population through a finite sample.
In theory, the only way to know the total species richness is to sample an infinite number of individuals. Without assumptions on the abundance distribution there can be an infinite number of arbitrarily rare species that a finite sampling experiment cannot possibly capture [18,5]. Recent work has expanded on this observation, helping to explain the apparent difficulty of constructing practical estimation methods by proving that no unbiased estimator exists in the general case [27,29]. Even non-parametric maximum likelihood estimators can suffer from extreme instability due to fitting abundance distributions with positive weight on arbitrarily small abundances [38,19,30].
The approach introduced by Harris [18] laid the groundwork for practical non-parametric moment-based estimation of species richness that others have built upon [7,8,9,28,29]. This general strategy relies upon a transformation of the problem of estimating the species richness to the problem of estimating an integral over an unknown distribution given its moments. These moments can be estimated as functions of the observed data and used to form a moment space. This moment space corresponds to the space of species abundance distributions that can explain the observed data. The associated integral has a maximum and minimum value over the constructed moment space, and these in turn correspond to maximum and minimum species richness that can explain the observed data. In agreement with previous results [29] we show that the maximum value is infinite regardless of the observed data. On the other hand the minimum value is finite and in theory will increase towards the true species richness as the dimension of the moment space increases, theoretically providing tighter lower bounds as more information is used.
The problem of estimating an integral over an unknown distribution given its moments is known as Gaussian quadrature, well studied in the numerical integration literature [11]. We use results from the theory of moment spaces to show that the solutions given by Gaussian quadrature are exactly the moment-based estimators of Harris [18]. The simplest case, when the moment space is degenerate, is precisely the Chao estimator [8]. Harris also presented a closed form solution when the moment space is two dimensional, but higher order solutions were previously unavailable. The use of Gaussian quadrature allows us to increase the dimension of the moment space and theoretically use an arbitrary amount of information to construct our estimators. We will show that this is complicated by the ill-conditioning of the procedure. In order to use the maximal amount of information while minimizing the effects of ill-conditioning, we develop an automated method to choose the order of the approximation. This alleviates some of the issues, but to fully avoid some of the severe instabilities we propose bootstrapping the observed histogram to obtain stable estimates.
The outline of the paper is as follows. We introduce the underlying model and the momentbased framework for estimating species richness in section 2. In section 3 we discuss how Gaussian quadrature can solve the moment-based problem. We follow with discussion and analysis of some of the methods for calculating Gaussian quadrature estimates and methods for combatting the inherent ill-conditioning of the problem in section 4. This includes an interesting counter-example specific to our problem to one of the prevailing suggested state of the art methods for numerical quadrature, as discussed in Laurie [23] and Gautschi [13]. In section 5 we compare the performance of the Chao estimator and our moment-based quadrature estimator on discrete abundance distributions. We illustrate the ill-conditioning problem exists even when the order of the approximation is chosen appropriately and how bootstrapping can mitigate the effects of ill-conditioning. Section 6 shows the performance of our estimator on simulated populations taken from large real world high-throughput sampling experiments. We close the paper with a discussion on directions for further investigation.

Moment-based estimation of species richness
Suppose that there are S species in the population. Let y i be the number of observations from species i, the so-called species count of i, for i = 1, . . . , S and let λ i be the corresponding Poisson rate. Assume that the rates are independently and identically distributed according to some distribution µ with support on the positive real line. The probability that j copies of a random species are sampled in the experiment is Note that conditional upon the total number of sampled individuals N = S i=1 y i , the species counts will follow a multinomial distribution with probabilities p i = λ i / S s=1 λ s . Define the count frequencies to be the number of species observed exactly j times in the sample By the compound Poisson assumption the expected value of n j is We assume no ordering of the species so that the observed data can be summarized by the count frequencies n 1 , n 2 , . . . and are a sufficient statistic regardless of S and µ.
The number of unobserved classes, n 0 , is unknown and has expectation Since S = N + n 0 , any estimate for n 0 immediately provides an estimate for S. Our objectve is simply to estimate n 0 , the unobserved species richness, using the observed data.

Estimating the unobserved species richness
Following ideas introduced by Harris [18], Chao [8] defined the measure ν such that dν(x) ∝ xe −x dµ(x). Note that this is a bijection with dµ(x) ∝ x −1 e x dν(x) and that both measures have the same support by definition. The moments of ν can be expressed as Since the observed count frequencies are always unbiased estimates of their expected values, we can substitute them in equation (1) to obtain estimated moments for ν. If we express E(n 0 ) in terms of dν, then Since the S term in E(n 0 ) has been absorbed into E(n 1 ) in equation (2) above, estimation of n 0 is free of S and equivalent to estimating the integral on the right hand side of the equation (2), The information we have on ν is expressed through the estimated moments. Therefore we can formulate our moment-based problem as follows: such thatν m = (m + 1)!n m+1 n 1 for m = 0, 1, . . . .
Estimates for the above integral, subject to the moment constraints, can be substituted in to equation (2) to obtain moment-based estimates of E(n 0 ). Since no unbiased estimator exists for the species richness in the general compound Poisson model [29], there is no estimator of ν that gives unbiased estimates of the integral (3). Therefore we will focus on upper and lower bounds, i.e. extremal estimates, of the integral subject to the moment constraints. These extremal estimates will correspond to lower and upper bounds for the species richness of populations that can explain the observed count frequencies, as reflected through the estimated moments.

Moment spaces and moment-based estimation
Consider the moment space of all distributions that satisfy the moment contraints The growth of the moments is bounded, i.e.ν m ≤ C(m + 1)! for some constant C. Therefore V M converges to a single point as M → ∞ [34, Proposition 1.5] and the measure ν is uniquely determined by its infinite sequence of moments. This in turn implies the identifiability of the species richness in the Poisson model, assuming one can accurately estimate all moments.
If dimension of the moment space is finite, i.e. M < ∞, as it will be in practice, then the space V M is convex, closed, and bounded [21]. The linear functional will therefore attain its extrema in the space V M . If V M is full rank, then the extrema will be acheived at unique discrete distributions of minimal degree contained in the moment space [18]. We can compute the extrema by finding the unique set of points 0 ≤ x 1 < · · · < x P and positive weights {w 1 , . . . , w P } that satisfy the moment constraints given by Lower extremal estimates can be obtained by taking M odd and P = (M + 1)/2. Upper extremal estimates can be obtained by taking M even, P = M/2, and fixing x 1 = 0. Given the points and weights that satisfy the system of equations (4), the estimated species richness is given by In the case of the upper extremal estimates, the estimated number of unobserved species is infinite, agreeing with the intuition that we can not exclude the possibility of an infinite number of species. If P = 1 and M = 1, then the system of equations is simply given by the two equations w 1 = 1 and x 1 =ν 1 = 2n 2 /n 1 . The corresponding estimate isn 0 = n 1 /x 1 = n 2 1 /(2n 2 ), exactly the estimator of Chao [8]. Similarly, Harris [18] derived an analytic solution to the system of equations when 2 points and 3 moments are used. Higher order analytic solutions are impossible as they involve finding the roots of degree 5 or higher polynomials.
If an estimator uses only the first few moments, and equivalently only the first few count frequencies, the estimator risks discarding a significant amount of the available information from the observed sampling experiment. In the case of highly complex populations, the underlying distributions may not be accurately explained by the first few moments. In modern data analysis the sampling experiments can be very large, allowing for accurate estimation of higher moments. We hypothesized that, conditional upon having accurate estimates for the higher moments, using more moments will lead to more accurate estimates of the species richness. This requires solving the system of equations (4). Since analytic solutions are impossible when higher moments are involved, we turn our attention to numerical solutions.
3 Numerical calculation of the extremal distributions with Gaussian quadrature The system of equations (4) is not unique to the problem of finding extremal values of an integral under moment constraints. We noted that the system is identical to the system of equations that Gaussian quadrature must satisfy [17]. Gaussian quadrature is a well studied numerical integration technique to approximate an integral over a unknown measure given its moments. This suggests an approach to use a larger number of moments by leveraging Gaussian quadrature techniques to solve the problem (3) numerically rather than analytically.

Gaussian quadrature
A numerical approximation, or quadrature, to the integral for some fixed integer P , and the approximation is exact for polynomials of degree 2P − 1 or less. The points x 1 , . . . , x P and weights w 1 , . . . , w P are collectively called the Gaussian quadrature rules.
Consider the standard basis polynomials: 1, x, . . . , x 2P −1 . Gaussian quadrature requires that This is exactly the system of equations (4) with M = 2P − 1. Note that any polynomial can be written as a finite sum of the standard basis polynomials. Therefore, if the system of equations (4) is satisfied then the quadrature rule will be exact for any polynomial of degree 2P − 1, and vice versa. Furthermore the Gaussian quadrature rules are unique [22]. This implies that the solution to the system of equations (4) is unique for fixed P , and the Gaussian quadrature rules and solutions to the system of equations given by Harris will correspond exactly. In the interpretation of Harris [18] and Mao & Lindsay [29] we are approximating ν with a discrete distribution based on the observed moments. Since the map from µ to ν is a bijection, we can interpret the quadrature rules as a moment-based discrete approximation to µ with increasing accuracy as the number of points increases (Figure 1 a). The corresponding estimates for the species richness will in theory converge from below to the true species richness (Figure 1

Choosing the order of the approximation
Choosing the order, or the number of points, in the quadrature estimator is similar to the problem of determining the order of a moment-based discrete approximation to a distribution. Of critical importance are the moment Hankel matrices, In previous applications [26,10], the observed data is typically considered to be independently and identically distributed according to the unknown distribution. In this case the moment Hankel matrices are constructed from the usual moment estimators, e.g. for i.i.d. observations x 1 , . . . , x n from an unknown distribution η,η m = 1 n n i=1 x i . This guarantees that the matrices will be nonnegative definite. The matrices are only singular when the support of the distribution is discrete with a support size smaller than the dimension of the Hankel matrices. One can then utilize the structure afforded by i.i.d. observations to construct statistical tests of non-negative definiteness and determine the order of the approximation [24,25].
In our problem, the moments are estimated from the count frequencies and not the usual moment estimators. This means that the moment Hankel matrix will not necessarily be positive semidefinite and we can not utilize the aforementioned tests. To ensure that the estimated quadrature rules correspond to a distribution on the positive real line, we require the estimated quadrature rules to be positive. This can be acheived by ensuring that all Hankel matrices up to order P + 1 are positive definite [29]. This requires the first 2P + 1 moments to be positive, while the quadrature rules only have 2P − 1 degress of freedom. To relax the requirement on the number of moments, we first define the shifted moment Hankel matrixH P = (ν i+j−1 ) 1≤i,j≤P . We next note that the quadrature rules will be strictly positive if and only if the generalized eigenvalues of the moment Hankel matrix and the shifted moment Hankel matrix are strictly positive [1]. In other words, all solutions z to the equation det H P − zH P = 0 are positive. A sufficient condition for this is when both H P andH P are positive definite [32]. To determine the number of points in the Gaussian quadrature estimator we test the determinants of H P andH P for successive values of P , stopping when either is not positive. Positive definiteness will be guaranteed when all lower order matrices have positive determinant by Sylvester's criterion. This requires the first 2P moments to be positive and allows us to use more information. In the worst case we can default to the single point case, exactly Chao's estimator.

Calculating Gaussian quadrature rules
To naively calculate the Gaussian quadrature rules one can use non-linear numerical solvers to calculate solutions to the system of equations (4) in a brute force manner. Numerical solvers have difficulty in converging in comparison to state of the art quadrature methods. These methods rely upon properties of the underlying distribution, particularly the associated orthogonal polynomials and their relationship with Gaussian quadrature rules [23].
Associated with a measure ν supported on the positive real line is a system of monic orthogonal polynomials {π i (x)} i≥0 . They are orthogonal in the sense that ∞ 0 π i (x)π j (x)dν(x) = 0 for i = j. These polynomials satisfy the three term recurrence for coefficients α i , β i > 0 and i = 0, 1, . . . Since π −1 (x) = 0, the value of β 0 is arbitrary but commonly set to 1. If the measure is uniquely determined by its moments then the three term recurrence is unique and completely characterizes the monic orthogonal polynomials [13]. We will show that the coefficients of the three term recurrence are intimately related to Gaussian quadrature rules so that estimation of the Gaussian quadrature rules can be broken into two parts: (1) estimating the coefficients of the three term recurrence for the moments and (2) estimating the quadrature rules from the estimated three term recurrence.
The mapping from the moments to the quadrature rules is known to be highly ill-conditioned [13]. Small changes in the moments can lead to extreme changes in the corresponding quadrature points and weights. This presents challenges for our approach, as variability in moment estimates can be substantial for higher moments. Most of the ill-conditioning arises from the first part above, while the second part is very well-conditioned [12].
Given the estimated coefficients of the three term recurrence we can construct the truncated Jacobi matrix, denoted J P . This is a symmetric tridiagonal matrix with {α 0 , α 1 , . . . , α P −1 } along the diagonal and { √ β 1 , . . . , β P −1 } on the main off-diagonals. The quadrature points x 1 , . . . , x P are the eigenvalues of J P and the quadrature weights w 1 , . . . , w P are the square of the first entry of eigenvectors. The eigenvalues and eigenvectors can be calculated simultaneously by the QR algorithm. Since J P is already tridiagonal, the usual Householder transformations are not required. Furthermore the full eigenvector is not needed. This allows for a modified QR algorithm that updates of the diagonal and off-diagonal elements of J P along with the first entry of the eigenvector at each iteration, without storing the full vectors, giving an algorithm with a time and space complexity linear in the number of moments used [14]. Since the truncated Jacobi matrix is symmetric, this procedure is perfectly conditioned.
Given the moments, the coefficients of the three term recurrence can be calculated by a number of algorithms. The first modern method was based on the Cholesky decomposition of the moment Hankel matrix [14]. The time complexity of this method is cubic in the number of moments. Due to exponential growth in the absolute condition number of the moment Hankel matrix [37], this method can suffer from ill-conditioning when using a large number of moments. To improve the conditioning a recursive method with quadratic time complexity was derived based on Chebyshev's work on orthogonal polynomials [36,33]. This method requires a known input measure, hopefully close to the unknown distribution, to modify the estimated moments. Accordingly this method is known as the modified Chebyshev algorithm. When no modifying measure is given we shall refer to it as the unmodified Chebyshev algorithm.
It is commonly suggested to use the modified Chebyshev algorithm with a known modifying measure and three term recurrence close to the unknown measure [2,23]. This is because if the modifying measure is serendipitously chosen to be exactly equal to the unknown measure then it will be recovered with no error. In previously considered applications the moments are assumed to be perfectly observed. The modified Chebyshev algorithm has not been investigated when the moments are not known exactly but instead are estimated with noise.
To evaluate the choice of the modifying measure we will establish a case where the coefficients of the three term recurrence can be derived analytically. We can then examine the situation where the true underlying measure is used as the modifying measure but only random estimates of the true moments are available, essentially if the modifying measure is chosen fortuitously.

The three term recurrence of negative binomial counts
Suppose that dµ is a Gamma density with shape parameter r and scale parameter p/(1 − p). In this case the number of observations from a random class follows a negative binomial distribution with size r and success probability p and dν(x) will be proportional to x r e −xφ with φ = 1/(1 − p).
We noted the similarity of the above measure to the measure that generates the classical monic associated Laguerre polynomials [35]. This motivates the following relationships to establish the orthogonal polynomials and the corresponding three term recurrence associated with the density dν. The proofs can be found in the Supplementary Materials.

Estimating the three term recurrence in the presence of error
To investigate the effect that error in estimating the moments has on the stability of the modified and unmodified Chebyshev algorithms, we simulated 1,000 independent samples from a population of one million classes, each class with counts that are independently distributed according to a negative binomial distribution with size parameter r = 1 and success probability p = 1/2. The three term recurrence was computed using the unmodified Chebyshev algorithm and the modified Chebyshev algorithm, with the modifying measure equal to the true underlying measure (dν ∝ xe −2x ) and the true three term recurrence as given in Proposition 2. The first few estimated coefficients of the three term recurrence (α 0 , α 1 , and β 1 ) were identical for both methods, indicating little effect of ill-conditioning in estimating the first few terms in the three term recurrence (Figure 2). For larger order coefficients, both methods exhibited increasing variance but the unmodified Chebyshev algorithm consistently exhibited less variance. We also observed a left shift in the distribution of the estimated coefficients away from the true values for the modified Chebyshev algorithm that was not present for the unmodified Chebyshev algorithm. This left shift may create problems downstream when estimating quadrature rules with a large number of points. In theory all of the coefficients of the three term recurrence are positive and we make this a strict requirement in our algorithm. The three term recurrence will be truncated and useful information will be thrown away unnecessarily, as exhibited in the lower percentage of positive quadrature rules in Figure 2 c.
The instability we observed in the modified Chebyshev algorithm can only be exacerbated when the true distribution is significantly different from the modifying measure [2]. The instability in calculating the coefficients of the three term recurrence will lead to instability in calculating the  quadrature rules. Small values of α and β can lead to smaller eigenvalues of the truncated Jacobi matrix, particularly for α since the trace of the truncated Jacobi matrix is equal to the sum of the eigenvalues. The small eigenvalues will then be inverted in the estimated integral, leading to large values in the estimated integral and instability in the estimated species richness (Figure 2). We found that in this case the performance of the unmodified Chebyshev algorithm was superior to the modified Chebyshev algorithm, even when the true measure is known, leading to more stable estimates of the species richness ( Figure 2 c and d).

Results
The performance of the Chao estimator (CE) has already been benchmarked on uniform populations and shown to perform well in comparison to other estimators [7, Chapter 6]. We examined the performance of the CE in the presence of heterogeneity and how using more information contained in the moments through our moment-based quadrature estimator (MQE) can improve the performance.

Two component discrete mixture
The simplest case of a perturbed uniform population is when there are two distinct uniform subpopulations. As long as both populations contribute to the first few count frequencies, the MQE should be able to distinguish the one point case from the two point. Examples where the twopoint might be indistiguishable from the one point is when few individuals are sampled from one subpopulation, one subpopulation is so prevalent that it contributes little to the first few count frequencies, or individuals in the two subpopulations have similar abundances.
For illustrative examples of two component discrete mixtures, free of the aforementioned problems, we consider populations of one million species with mixtures of 40%, 30%, 20%, and 10% of the population represented five and ten times more abundant than the remaining species, described in Table 1. We refer to the rarer species as the background and the more abundant species as the foreground. We sampled an average of one million individuals from each population one hundred times and calculated the CE and MQE (Figure 3).
In all samples the MQE correctly selected a two point approximation and the estimates were centered around the true species richness. The CE's were all lower bounds, as they should be in the presence of heterogeneity [3]. The average improvement of the MQE to the CE ranged from 344,250, or about one third of the population, to 38,443, less than five percent of the population. The difference between the two estimators tended to be less when the foreground population was smaller. In large part, this is because the majority of the foreground population in these cases is observed and contributes little to smaller count frequencies. For example, for population h, shown in Figure 3h, 97% of the foreground population is expected to be observed. It contributes approximately 3% of the singletons and 9% of the doubletons in expectation but about half of the species in n 3 . Thus the first two count frequencies closely reflect the fact that the vast majority of the unobserved species are uniform and the CE will represent this. In contrast, in population b about half the singletons are expected to be from the foreground. The CE has a hard time in estimating the unobserved species due to the distortion from the foreground in the first two count frequencies (Figure 3b). One of the main advantages to the use of CE in estimating species richness is its stability. Even in the ideal equiprobable case the CE performs well compared to maximum likelihood estimators [7,Chapter 6]. The addition of extra information in estimating the species richness through n 3 and n 4 introduces extra variability into the MQE compared to the CE. Obviously there is a bias-variance tradeoff, where estimation using only the first two count frequencies, in the CE, is more biased but less variable while estimation using the first four count frequencies, in the MQE, is less biased but more variable. When the results are measured by root mean square error (RMSE) the MQE shows an order of magnitude lower RMSE than the CE.

Higher order mixtures: three and four components
The next obvious step is to examine higher order mixtures, specifically three and four component mixtures. We investigated four populations, each composed of one million species. The first two populations had rate distributions that followed a three component mixture, The improvement of the MQE over the lower order estimates naturally brings more variance into the estimates due to the use of more information through more moments. The ill-conditioned nature of the estimation compounds this difficulty, so that some of the MQEs are much, much larger than the true species richness. The effect of this is to typically decrease the performance of the higher order approximations in terms of the root mean square error, though the tradeoff between the number of moments used in the estimates and the RMSE depends on the properties of the population.
To combat the ill-conditioning issue we suggest bootstrapping the observed counts frequencies, which conveniently follow a multinomial distribution [39]. Averaging over the bootstrapped MQE (bMQE) reduces the effect of those estimates that suffer from ill-conditioning and significantly improves the performance of the moment-based quadrature estimator. This results in the bMQE having the highest average species richness estimates and lowest RMSE (Figure 4e-h).

Examples
It is always the desire of researchers to test methods on real populations. A few researchers, such as Carothers [6] have designed experiments where the properties of the population are completely a b c d  known. Unfortunately such cases are few and small in scale by necessity. Our focus is on large scale sampling experiments where the higher order moments are estimated with more accuracy and contain information that can significantly improve the estimator. In order to simulate large heterogenous populations with known properties and are similar to populations encountered in real world sampling experiments we took a numerical simulation approach. We took large observed experiments as the population and sampled with replacement from the experiment to simulate sampling individuals. We examined several situations where species richness estimators have been applied and used for inference, specifically high-throughput RNA-seq, immune repertoire, and metagenomics.
To simulate high-throughput RNA-seq experiments we took an extremely deep RNA-seq experiment of 1,095,239,309 mapped single end 100 base pair reads from a murine model of acute myeloid leukemia. We obtained 74,722,599 distinct reads, 173,485 observed RefSeq coding exons, and 402,196 unique junctions. The total number of observed exons and junctions was resepctively 820,653,733 and 235,107,245. The coefficients of variation are 39.8, 5.1, and 8.9 for the read (Lib-Size), exonic (Exons), and slice junction (Junctions) populations. For the immune repertoire we examined the combined T-cell receptor β repertoire of 39 individuals of varying ages [4], using the inferred amino acid sequence to determine the repertoire (TCR). The combined data contained 16,704,121 TCRβ species with 38,701,991 total observations and a CV of 66.8. To simulate metagenomic sampling experiments we took α-diversity source data from the MG-RAST database [31]. We examined two populations: soil microbiome of the Penang Mangrove forest (PenangMangrove, experiment ID 4510219.3) and the combined population of 14 marine sediment microbiome samples from the Gulf of Mexico following the Deepwater Horizon oil spill (GulfSediment, project ID 3023) These contained, respectively, 69,454 and 41,720 annotated species with 247,106,874 and 577,099,067 total observations and CVs of 11 and 11.1. All of the observed count frequencies are available in the Supplementary Materials.
For each population we examined three levels of sampling, varying across three orders of magnitude. We selected these levels based on existing experiments to approximately correspond to shallow, moderate, and deep sampling. For the RNA-seq experiments we sampled levels equivalent to 500,000 (500K), 5,000,000 (5M), and 50M mapped reads. In the metagenomic sampling experiments we sampled 10K, 100K, and 1M annotated species. Finally, in the T-Cell repertoire we sampled 100K, 1M, and 10M T-Cell β sequences (TCS). At each level we took 100 independent samples from the constructed samples.
We measured the average performance of the Chao estimator and the moment-based quadrature estimators with the median species richness estimate (S). We also considered the variability of the estimates, quantified by the standard error (SE). There is a tradeoff to consider between the average estimate and variability. For example, one can use the observed number of species as an estimate of species richness. This simple estimator will have low variability but is typically a very poor estimator. To capture the bias-variance trade-off we used the root mean square error (RMSE). These results are summarized in Table 2.
The moment-based quadrature estimators were larger than the Chao estimator in all samples and in only 3 out of the 600 deepest samples was the MQE larger than the true species richness. The CE, M2QE, and M3QE always underestimated the true species richness, as did the MQE at the lower sampling depths. This method of simulation essentially constructs high, but finite, dimensional mixtures. In these cases, using just a few moments is not sufficient to estimate the species richness. The use of more moments, at least through our framework, can only increase the estimated species richness. But this naturally increases the variance and we observe increasing variability in the maximum number of moments used. The question that naturally follows is if the increase in the estimated species richness is worth the larger variability.
In all cases the MQE had the lowest RMSE, indicating that the increasing accuracy of using more points is worth the increase in variance. Sometimes the improvement was small, as was the case in the metagenomic samples, while in some cases the MQE was an order of magnitude larger than the CE, as was the case in the Exons population with a 500K read sample. Interestingly, for the metagenomic populations (GulfSediment & PenangMangrove) the MQE and M3QE were exactly equal in all but a few cases. This is most likely due to the fact that the MQE used more than 3 points extremely infrequently for the estimates in this case.
We noted the estimated species richness for the metagenomic populations is several times lower than the true species richness at the deepest level of sampling. Even when sampling at an order of magnitude larger (10M sequences from annotated species), the MQE still significantly underestimates the metagenomic species richness. This may indicate that larger experiments are needed to accurately quantify the rarest species in metagenomic sequencing experiments. On the other hand, it is clear that deep sampling or sequencing is required to capture transcriptional diversity, including both exonic and transcriptional diversity, and T-Cell repertoire diversity. It may be that our analysis ignores extremely rare sequences missing from the original population so that even deeper sampling is required in practice.

Discussion
The general moment-based framework for estimating species richness was developed over 50 years ago by Harris [18] and later explored by Chao [7], yet the practical application has been impeded two primary issues: the lack of deep sampling experiments where higher moments can be accurately estimated, and efficient algorithms to compute the estimates. Inspired by the recent explosion of deep sampling experiments brought on by technological improvements, we presented algorithms specifically designed to calculate species richness in the moment-based framework of Harris to use the maximum amount of information possible from the observed species counts. We showed that our method improves upon the Chao estimator in the presence of heterogeneity and will default to the Chao estimator in the worst case.
There have been a multitude of new approaches in recent years for estimating species richness, including penalized non-parametric maximum likelihood estimation [39] and non-parametric non-linear regression of the count frequency ratios [40]. The moment-based quadrature estimator provides a new avenue of research into estimating species richness. Possible directions of future research include investigating algorithmic improvements for estimating the three term recurrence in the presence of error, a woefully under explored area of research that we briefly discussed in this paper. Such work to combat the ill-conditioning in the algorithm and improve the stability of the estimates will facilitate the use of more moments in the estimator. Our work suggests new conditions for smoothing count frequencies, a topic with a long history in the species richness problem [15,16]. Rather than smoothing the count frequencies directly, one can smooth the moments through the moment Hankel matrices to ensure that the moments correspond to a true distribution, similar to the problem of finding the nearest covariance matrix [20], but with the added condition of constant off-diagonals.
The development of new sampling technology, particularly in high-throughput sequencing, is increasing our understanding of large and highly heterogeneous populations. The development of fast, efficient, and accurate algorithms will help further our understanding of these populations. The method we presented represents a significant improvement on previous non-parametric moment based estimators and allows for fast and accurate estimates of species richness. We have made the software available as part of the open source preseq package (http://smithlabresearch.org/ software/preseq/).
We will show that these polynomials are the only system of polynomials orthogonal under the measure dν(x) = x r e −xφ for r > 0 and φ > 0. To this end first note that the polynomial can be rewritten in the Rodrigues' formula Consider the classical associated Laguerre polynomials [35], orthogonal under the measure dω(x) = x r e −x dx and given by By othogonality the following holds for all positive integers i = j, x j+r e −x x r e −x dx = 0.
Consider the substitution x = φy with dx = φdy. Then the left hand side of equation (11) above is equal to By equation (11), the above integral in equation (12)  Note that trivially ∞ 0 L (r,φ) i (y) 2 dν(y) > 0 for all non-negative integers n. Therefore the system of monic polynomials {L (r,φ) i (x)} i≥0 is orthogonal under the measure dν(x) = x r e −xφ dx. Since the monic polynomials can be constructed by the Gram-Schmidt algorithm, uniqueness follows immediately from the positive definiteness of the measure dν.

A.2 Proof of Proposition 2
Proof. Since {L (r,φ) i (x)} i≥0 is a system of monic polynomials orthogonal under the measure dν, {L (r,φ) i (x)} i≥0 must satisfy the three term recurrence of equation (7) for some α i > 0 and β i > 0. Therefore, the following holds for all i ≥ 0, Substituting L (r,φ) i (x) defined in equation (8) results in the following equation Since equation (14) must hold for all x, the reduced coefficients of both sides must be equal. We can group terms on both sides of the equation according to the degree of x to solve for α i and β i . These polynomials are monic, so the highest order term is x i+1 . The term of second highest order, as written on both sides of the equation, simplifies to Γ(i + r + 2)(i + 1)! Γ(i + 1 + r)i! φ −1 (−1)x i = Γ(i + r + 1)i! Γ(i + r)(i − 1)! φ −1 (−1) −1 x i − α i x i .
Further simplifying the coefficients of x i results in (i + r + 1)(i + 1) = (i + r)i + α i φ, providing α i = (2i + r + 1)/φ. Now consider the lowest order terms. The constants written on both sides of equation (14) are After simplification, and substitution with α i above, we arrive at β i = (i 2 + ir)/φ 2 .
The above calculated α i and β i , obtained using the second highest and lowest order coefficients, can be used to verify the equality of equation (14) for all other coefficients.