The distribution of the number of mutations in the genealogy of a sample from a single population

The number $K$ of mutations in the genealogy of a sample of $n$ sequences from a single population is one essential summary statistic in molecular population genetics and is equal to the number of segregating sites in the sample under the infinite-sites model. Although its expectation and variance are the most widely utilized properties, its sampling formula (i.e., probability distribution) is the foundation for all explorations related to K. Despite existence of an analytic sampling formula, its numerical application is limited due to susceptibility to error propagation. This paper presents a new sampling formula for $K$ in a random sample of DNA sequences from a neutral locus without recombination, taken from a single population evolving according to the Wright–Fisher model with a constant effective population size, or the constant-in-state model, which allows the effective population size to vary across different coalescent states. The new sampling formula is expressed as the sum of the probabilities of the various ways mutations can manifest in the sample genealogy and achieves simplicity by partitioning mutations into hypothetical atomic clusters that cannot be further divided. Under the Wright–Fisher model with a constant effective population size, the new sampling formula is closely analogous to the celebrated Ewens’ sampling formula for the number of distinct alleles in a sample. Numerical computation using the new sampling formula is accurate and is limited only by the burden of enumerating a large number of partitions of a large K. However, significant improvement in efficiency can be achieved by prioritizing the enumeration of partitions with a large number of parts.