Simulating probability distributions in the coalescent (Q1340495)

We describe some computational algorithms for computing probability distributions for sample configurations from the finite-sites models in population genetics. Our main interest is the development of computational methods for estimating substitution rates for DNA sequence data using maximum likelihood. In broad terms the method we propose has two main aspects: The first is the derivation of a recursion satisfied by the sampling distribution of the data -- the distribution of the number of each allele observed in the sample. This uses the structure of coalescent that describes the ancestral relationships among the individuals in the sample. The recursion is, typically, too large a system to solve by direct numerical methods. We adopt instead a simulation approach. This leads to the second aspect, the construction of a Markov chain with a set of absorbing states in such a way that the required sampling distribution is the mean of a functional of the process up to the absorption time. This provides a conceptually simple framework for simulating the likelihood of the data for given parameter values, or, more generally, for a set of parameter values. The method is particularly attractive in practice: it is simple to program, and can be extended to cover other features of interest such as the infinitely-many-sites process, recombination, selection, and variable population size.