ON A BAYESIAN APPROACH TO ESTIMATING CLASS NUMBER

STEPHEN MONTGOMERY-SMITH

Abstract. We examine a Bayesian approach to estimating the number of classes in a population, in the situation that we are able to take many independent samples from an infinite population.

Introduction. Suppose we capture 1000 birds from an effectively infinite population, and observe 15 different species. Estimate how many species have not been observed. While it is clear that one can only get a lower estimate (because there may be a very large number of extremely rare species which we have essentially no hope of observing), nevertheless if, say, 4 of the species appear only once, we would assume that most likely that there are species not observed, because if 15 were the exact number we would likely have not observed every one of them.

More generally, suppose that we have an observation O of n independent samples from an infinite population. Each sample belongs to one of M classes C_{1}, C_{2}, …, C_{M}, and the probability that a sample belongs to class C_{i} is p_{i}, where 0 ≤ p_{i} ≤ 1 and ∑_{i=1}^{M} p_{i} = 1. The goal of this paper is to present a Bayesian method to estimate the class number, M. The difficulty is to estimate the number of classes that have not been observed.

A survey of various approaches is given in [2]. In particular, we point out a rather effective method due to A. Chao [1]. If m_{k} is the number of classes observed k times, then a lower estimator for the number of classes not observed is m_{1}^{2}/2m_{2}.

The Bayesian approach will be to treat M and p = (p_{1}, …, p_{M}) as random variables with prior distributions. After computing the posteriors, we will integrate out p. This has already been explored, for example, in [3], where the final distribution is obtained using Monte-Carlo techniques. In this paper we present some prior distributions for which we are able to give explicit formulae for the final distributions.

The purpose of this paper is not to advocate for the Bayesian approach, but rather to identify that in this particular situation that the Bayesian approach seems to give rather poor quality answers. We will present a computer generated simulation in which our methods will significantly overestimate. In particular, the beta prior distributions will be shown to give significantly bland answers, because the final formulae do not depend upon the numbers m_{k}, but only upon m = ∑_{k=1}^{n} m_{k}, the total number of classes observed, and n = ∑_{k=1}^{n} km_{k}, the sample size.

The posterior. Let n_{1}, n_{2}, …, n_{m} denote the number of samples from each observed class, that is, a sequence in which each number k appears exactly m_{k} times.

Now we do not know which class is which, and so

Pr(O|(M, p)) = ∑_{σ∈B_{m,M}} ∏_{i=1}^{m} p_{σ(i)}^{n_{i}},

where B_{m,M} denotes the set of one to one mappings from {1, 2, …, m} to {1, 2, …, M}. (We mention that this formula assumes that we are also told the order in which the samples are picked. If we do not know the order in which they are picked, and we are simply given the numbers n_{1}, n_{2}, …, n_{m}, then there is an extra factor n!/n_{1}! ⋯ n_{m}!m_{1}! ⋯ m_{n}!. However this factor is independent of M and so it will disappear in the final normalization that takes place in Bayes’ Theorem.)

The prior for p given M is a distribution μ_{M} on the simplex

Δ_{M−1} = {p = (p_{1}, …, p_{M}) : p_{i} ≥ 0, ∑_{i=1}^{M} p_{i} = 1} ,

which, since we do not know a priori which class is which, must be invariant under rearrangements of p_{1}, …, p_{M}. Thus

Pr(O|M) = ∫_{Δ_{M−1}} Pr(O|(M, p) dμ_{M}(p) = \frac{M!}{(M − m)!} ∫_{Δ_{M−1}} ∏_{i=1}^{m} p_{i}^{n_{i}} dμ_{M}(p).

Prior distributions from populations. The prior distributions of p that we will consider will be derived from the size of the population of each class, X_{i}, which are very large (that is, effectively infinite) random variables, but independent and identically distributed. We will assume that each X_{i} has a continuous probability distribution λϕ(λx) dx, where λ > 0 is very small. In the appendix we will show that

Pr(O|M) = \frac{M!}{(M − m)!(n − 1)!} ∫_{s=0}^{∞} s^{−M−1}h_{0}(s)^{M−m}h_{1}(s)^{m_{1}}h_{2}(s)^{m_{2}} ⋯ h_{n}(s)^{m_{n}} ds,

where

h_{k}(s) = ∫_{0}^{∞} x^{k}e^{−x}ϕ(x/s) dx.

Beta prior distributions. This distribution arises when the population size of each class is given by the Gamma distribution, namely the distribution is ϕ(x) = x^{a}e^{−x}/Γ(a + 1), where a > −1. Then dμ_{M}(p) is a beta distribution, that is, it is proportional to (p_{1} ⋯ p_{M})^{a} times Lebesgue measure on Δ_{M−1}. The case a = 0 is the uniform distribution on Δ_{M−1}. A direct calculation shows

∫_{Δ_{M−1}} ∏_{i=1}^{M} p_{i}^{n_{i}} dμ_{M}(p) = \frac{Γ(M + Ma + 1)}{Γ(M + Ma + n + 1)} ∏_{i=1}^{m} \frac{Γ(n_{i} + a + 1)}{Γ(a + 1)}.

After applying Bayes’ Theorem, and remembering the normalization that takes place, we find that

Pr(M|O) ∝ \frac{M!Γ(M + Ma + 1)}{(M − m)!Γ(M + n + Ma + 1)} Pr(M),

where the constant of proportionality is so that the probabilities on the left hand side add up to 1. We notice that the answer does not depend upon the values of n_{1}, n_{2}, …, n_{m}, but only upon the sample size n, and the number of classes observed m.

Uniform distribution on the sphere. Here we suppose that the population sizes of each class have a semi-Gaussian distribution, that is ϕ(x) = \sqrt{2/π} e^{−x^{2}/2}. This gives rise to a distribution which is formed by projecting the normalized uniform distribution on the positive sphere

S_{M−1}^{+} = {q = (q_{1}, q_{2}, …, q_{M}) : q_{i} ≥ 0, ∑_{i=1}^{M} q_{i}^{2} = 1}

onto Δ_{M−1} via the map q ↦ q/(q_{1} + q_{2} + ⋯ + q_{M}). In this case we obtain

h_{k}(s) = \sqrt{\frac{2}{π}} ∫_{0}^{∞} x^{k}e^{−x} exp(−x^{2}/2s^{2}) dx.

Although h_{k}(s) can be computed explicitly as (−1)^{k}s^{k+1} times the kth derivative of exp(s^{2}/2) erfc(s/\sqrt{2}), the resulting formula is numerically unstable for large s. Instead it is rather easily computed by a standard numerical integration package.

Uniform distribution on the rectangle. This next possibility is if the population of each class is uniformly distributed. In this case h_{k}(s) is a lower incomplete Gamma function, which can easily be numerically computed.

A numerical simulation. A computer simulated a sample of n = 500 taken from a population of 100 classes, of which 50 of the classes occurred with a probability of 0.05, and the other 50 classes occurred with a probability of 0.15. The data obtained had m_{1} = 12, m_{2} = 15, m_{3} = 12, m_{4} = 8, m_{5} = 6, m_{6} = 15, m_{7} = 8, m_{8} = 4, m_{9} = 4, m_{10} = 8, m_{11} = 2, m_{14} = 1, m_{15} = 2, and all other m_{k} = 0. Thus the number of observed classes is m = 97. The Chao estimate gives 3.2 for the number of unobserved classes, which is very close to the correct answer.

Using the improper prior distribution where Pr(M) is constant, then the beta with a = 0, uniform spherical, and uniform rectangular priors for p respectively give a final probability of M ≤ 100 equal to 3.46 × 10^{−6}, 0.00028, and 0.0048, and a mid-90th percentile range for M of 121–131, 106–121, and 103–115. It can be seen that the uniform rectangular prior seems to work the best, but all of them significantly overestimate the class number.

Appendix: derivation of prior distribution from populations. We have that

∫_{Δ_{M−1}} ∏_{i=1}^{M} p_{i}^{n_{i}} dμ_{M}(p) = ∫_{ℝ_{+}^{M}} \frac{x_{1}^{n_{1}}x_{2}^{n_{2}} ⋯ x_{m}^{n_{m}}}{(x_{1} + ⋯ + x_{M})^{n}} λ^{M}ϕ(λx_{1}) ⋯ ϕ(λx_{M}) dx.

Since

∫_{s=0}^{∞} s^{n−1}e^{−sy} ds = \frac{(n − 1)!}{y^{n}},

the right hand side may be seen to equal

\frac{λ^{M+n}}{(n − 1)!} ∫_{s=0}^{∞} s^{n−1} ∫_{ℝ_{+}^{M}} x_{1}^{n_{1}}x_{2}^{n_{2}} ⋯ x_{m}^{n_{m}}e^{−λ(x_{1}+···+x_{M})s}ϕ(λx_{1}) ⋯ ϕ(λx_{M}) dx ds,

which after making the substitution x_{i} ↦ x_{i}/λs can be written as

\frac{1}{(n − 1)!} ∫_{s=0}^{∞} s^{−M−1}h_{n_{1}}(s) ⋯ h_{n_{m}}(s)h_{0}(s)^{M−m} ds.

Acknowledgments. The author would like to extend his sincere thanks to Anthony Garrett introducing him to the Bayesian approach to this particular problem. He also extends his gratitude to Frank Schmidt for bringing the class number problem and reference [1] to his attention.

References

[1] Chao, A. 1984. Non-parametric estimation of the number of classes in a population. Scandinavian Journal of Statistics 11, 265-270.

[2] Chao, A. 2005. Species richness estimation, Pages 7909-7916 in N. Balakrishnan, C. B. Read, and B. Vidakovic, eds. Encyclopedia of Statistical Sciences. New York, Wiley.

[3] Rodrigues J., Milan L.A., Leite J.G. 2001. Hierarchical Bayesian estimation for the number of species. Biometrical J., 43, 737-746.

Math Dept, University of Missouri, Columbia, MO 65211, U.S.A.

  E-mail address: stephen@math.missouri.edu