# A matrix-analytical sampling formula for time-homogeneous coalescent processes under the infinite sites mutation model

Asger Hobolth, Simon Boitard, Andreas Futschik, Raphael Leblois

► To cite this version:

Asger Hobolth, Simon Boitard, Andreas Futschik, Raphael Leblois. A matrix-analytical sampling formula for time-homogeneous coalescent processes under the infinite sites mutation model. Theoretical Population Biology, 2025, 163, pp.62-79. 10.1016/j.tpb.2025.03.002. hal-05030389v2

HAL Id: hal-05030389

https://hal.inrae.fr/hal-05030389v2

Submitted on 25 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# A matrix-analytical sampling formula for time-homogeneous coalescent processes under the infinite sites mutation model

Asger Hobolth<sup>a, ∗</sup>, Simon Boitard<sup>b</sup>, Andreas Futschik<sup>c</sup>, Raphael Leblois<sup>b</sup>

ARTICLE INFO

Keywords:

Coalescent theory

Infinite sites mutation model

Phase-type theory

Poisson thinning

Probability Generating Function

Sampling formula

A B S T R A C T

In this paper we develop a general framework for calculating the probability of a genetic sample under a time-homogeneous coalescent process and the infinite sites mutation model. The evolutionary model that we consider can be characterized as a two-step procedure: A coalescent process that describes the ancestral relatedness of the samples and a sprinkling of mutations in separate sites on the ancestral tree according to a Poisson process. The coalescent process is defined using multivariate phase-type theory. The requirements are a rate matrix that determines the transition rates between the ancestral states, an initial state probability vector, and a reward matrix that informs about the characteristics of the ancestral states. For example, the reward matrix could contain information about the number of singleton, doubleton or higher-order lineages in the ancestral states. We analyze the probability generating function for the evolutionary model as a function of the initial state probability vector, the transition rate matrix, the reward matrix, and the mutation rate. The matrix-analytical expression of the probability generating function allows us to develop a general method for calculating the probability of a population genetic data set. We demonstrate that the method is computationally attractive for a small number of mutations and provide a simple and easy-to-implement algorithm for determining the probability of a sample from the evolutionary model. The method is computationally stable and only involves a single inverse matrix operation, matrix multiplications and matrix additions. We provide comprehensive understanding of the procedure by detailed calculations and discussions of several elementary examples. These examples include different sample representations (labeled samples and the site frequency spectrum) and different demographic and genetic models (the structured coalescent and the Beta-coalescent). We apply the sampling formula to calculate probabilities of spectra for the Kingman coalescent and the Beta-coalescent. Even for a small number of samples and mutations we find that the probabilities for spectra vary in huge orders of magnitudes. We compare the probabilities of the spectra to the values of Tajima's *D*-statistics, and find that the *D*-statistic is a poor predictor for the probability of a spectrum. Finally, we investigate how the probabilities of the spectra vary with the parametrization of the Beta-coalescent.

# 1. Introduction

The purpose of this paper is to develop a general framework for calculating the probability of a genetic sample under a time-homogeneous coalescent process and the infinite sites mutation model. In the infinite sites model, mutations are sprinkled on branches according to an independent Poisson process. We have previously shown (Hobolth et al., 2019) that the branch lengths in the ancestral process of a time-homogeneous coalescent model are distributed according to Kulkarni multivariate phase-type distribution (Kulkarni, 1989). The joint distribution function of the branch lengths in Kulkarni multivariate phase-type distribution is not available in closed analytical form (Bladt and Nielsen, 2017). However, fortunately, the Laplace transform for the multivariate phase-type distribution is available in a rather attractive matrix-analytical form. We apply this property in the derivation of our matrix-analytical sampling formula for the infinite sites model. In particular we take advantage of the closed matrix-analytical form of the probability generating function (PGF): we demonstrate how to differentiate the PGF to obtain the probability of a sample

<footnote><sup>a</sup> Department of Mathematics, Aarhus University, Aarhus, Denmark</footnote>

<footnote><sup>b</sup> CBGP, INRAE, CIRAD, IRD, Institut Agro, Univ Montpellier, Montpellier, France</footnote>

<footnote><sup>c</sup> Institute of Applied Statistics, Johannes Kepler University, Linz, Austria</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: asger@math.au.dk (A. Hobolth), simon.boitard@inra.fr (S. Boitard), andreas.futschik@jku.at (A. Futschik), raphael.leblois@inrae.fr (R. Leblois).</footnote>
from a time-homogeneous coalescent process under the infinite sites mutation model. We argue that our sampling formula has a nice probabilistic interpretation in terms of thinning the transition probability matrix for the total number of mutations.

The probability generating function for a general multivariate phase-type distribution with independent Poisson distributions was derived using the phase-type framework in Hobolth et al. (2021), where it was used to determine the distribution of linear functions of the site frequency spectrum. Previously, the probability generating function (PGF) has been applied extensively in the work by Robert Griffiths and Simon Tavare going back to Griffiths (1982). Much of their work was recently summarized in Griffiths and Tavare (2018).

Konrad Lohse and Nick Barton have also applied the PGF procedure in a series of papers on likelihood inference in the structured coalescent. Their first work is Lohse et al. (2011), and they are still developing the methodology. Similar to our approach they are taking suitable partial derivatives of the PGF and evaluating them at zero in order to calculate the probability of a sample. Lohse et al. (2011) compute the derivatives symbolically using Mathematica (2024). This can also be done for different demographic models such as the isolation with migration model. Although some symmetries can be exploited, symbolic computations are only possible for small samples of up to six sequences, and a small number of mutations. To keep the number of mutations small, the authors propose to compute block-wise likelihoods for sufficiently short DNA segments.

Recently, Bisschop (2022) proposed a more efficient way of computing the generating function using a graph traversal algorithm that relies on the state space graph. This simplifies the derivation of the joint distribution of branch lengths. However, computing higher-order derivatives of the PGF seems to remain a challenge also with this approach unless the number of mutations is small.

Marcy Uyenoyama and co-authors have also extensively explored the PGF framework. Uyenoyama et al. (2019) derived the PGF for the allele frequency spectrum of the structured coalescent. The framework was based on earlier work on generating functions for structured populations and a panmictic population, and in particular (Uyenoyama and Takebayashi, 2004; Ganapathy and Uyenoyama, 2009; Kumagai and Uyenoyama, 2015). These papers also contain several references to earlier work on the PGF for the basic coalescent process with the infinite alleles or infinite sites mutation models.

One possible outcome of the PGF approach is to derive the probability distribution of the site frequency spectrum (SFS) at any given non-recombining locus. The SFS has a long history in population genetics, and has been used to study coalescent models with variability in population size (e.g. Griffiths and Tavare, 1998; Bhaskar et al., 2015; Liu and Fu, 2015), multiple mergers (e.g. Spence et al., 2016), and population structure (e.g. Chen, 2012; Gutenkunst et al., 2009; Excoffier et al., 2013). While several of these approaches rely on analytical derivations under the coalescent and on state spaces that are similar to the ones considered here, note that they all focus on the expected SFS obtained from mutations at independent loci. This is in stark contrast to our approach as well as the other PGF approaches mentioned above. In particular, the PGF approach carries richer information about the topology at each specific locus, which should improve parameter inference.

Our procedure is formulated on the matrix-analytical level, and we use reward-transformations (Kulkarni, 1989) of the ancestral process to capture the dependencies of branch lengths in the history of a sample (Hobolth et al., 2019). We find that working on the level of matrices allows for easy implementations and using reward-transformations substantially increases a detailed understanding of the formulas. In particular we provide a natural interpretation in terms of Poisson thinning of mutations.

In general our framework takes advantage of the close connection between queuing theory and coalescent theory (Hobolth et al., 2024). This connection dates back to at least Kingman (1982). The celebrated Ewens Sampling Formula (Ewens, 1972) is now more than fifty years old. A fascinating account of the history and importance of the formula is given in Rosenberg (2020). Ewens Sampling Formula gives the probability of a sample of alleles. John Kingman worked on queuing theory in the 1960's, and Geoff Watterson suggested a possible relation between the infinite alleles mutation model and queuing theory. This led John Kingman to develop the *n*-coalescent (Kingman, 1982). Section 3.5 in Tavaré (2004) derives Ewens Sampling Formula from the *n*-coalescent and under the infinite alleles mutation model. It is perhaps not surprising that a framework originally developed for queuing theory, namely Kulkarnis multivariate phase-type distribution (Kulkarni, 1989), can also be applied to derive a sampling formula for time-homogeneous coalescent processes under the infinite sites mutation model.

Our paper is organized as follows. In Section 2 we provide the theoretical background for our matrix-analytical sampling formula for the infinite sites model. It is not an easy task to get familiar and comfortable with the quantities that enter into the PGF. Furthermore, the 'Laplacian curtain' (Neuts, 1981) makes it hard to develop an intuitive understanding of the formulas. We try to remedy this situation by supplementing the general theory with an example of the site frequency spectrum with sample size three. Furthermore, in this section, we ask for persistence and patience from the reader. In Section 3 we provide three more illustrative examples and elaborate on the example from Section 2. The general sampling formula for the probability of a sample is provided in Section 4.1. The sampling formula is constructive in the sense that it immediately suggests an accompanying algorithm, which we describe in Section 4.2. In Section 5 we provide additional elementary examples: a labeled sample of size three and a labeled two-island model of sample size two. Section 6 is of a more applied nature. In particular, we investigate the sampling probabilities of site frequency spectra with a fixed number of mutations and a fixed sample size. We ask: is Tajima's *D*-statistic (Tajima, 1989) a good predictor of the sampling probability? And we ask: how do the probabilities of the site frequency spectra change with the parametrization of the Beta-coalescent? The implementation of the sampling formula and the code used to generate the figures in the paper are available at https://github.com/AsgerHobolth/MatrixAnalyticalSamplingFormula. The paper ends with a discussion about improved computational efficiency, relaxation of the homogeneity assumption of the coalescent process, and the prospects of improved Monte Carlo-based inference in population genetics based on the new sampling formula.

## 2. Probability generating function for the infinite sites model

In this section we provide the theoretical background for our matrix-analytical sampling formula for the coalescent process under the infinite sites mutation model. In particular we describe the probability generating function (PGF) for a time-homogeneous coalescent process with Poisson mutations in separate sites. In order to calculate the probability of a sample configuration we need to differentiate the PGF an appropriate number of times with respect to each coordinate and evaluate at zero. The contribution of this paper is to provide a general procedure for accomplishing this task.

Consider the joint distribution of *d* branch lengths ($B_1, \dots, B_d$) in a time-homogeneous coalescent model (see e.g. Wakeley, 2008, for an introduction to coalescent theory). For example, the branch lengths of interest could be the stochastic vector ($B_1, B_2, B_3$) with entries corresponding to the total singleton branch length $B_1$, total doubleton branch length $B_2$, and total branch length subtending three leaves (total tripleton branch length) $B_3$, in a standard coalescent model with a sample size of *n* = 4. A realization from the coalescent process results in a genealogical tree. Associated with a given tree is a collection of time intervals and the number of leaves subtended by extant lineages within each time interval.
Initialization $\alpha = (1,0)$ and sub-intensity matrix $S = \begin{pmatrix} -3 & 3 \\ 0 & -1 \end{pmatrix}$

Time to absorption (most recent common ancestor) $\tau = T_3 + T_2 \sim \text{PH}(\alpha, S)$

State of the coalescent process before absorption $\{X_t\}_{0\le t\le \tau}$

Reward transformation: singleton branch length $B_1 = 3T_3 + T_2 = \int_0^\tau r_1(X_t)dt$, where $r_1 = (3,1)$

Reward transformation: doubleton branch length $B_2 = T_2 = \int_0^\tau r_2(X_t)dt$, where $r_2 = (0,1)$

<figure><img src="image_4.png" /></figure>

Joint distribution $B = (B_1, B_2) \sim \text{MPH}(\alpha, S, R)$, where reward matrix $R = (r_1, r_2) = \begin{pmatrix} 3 & 0 \\ 1 & 1 \end{pmatrix}$

Total branch length $B_{\text{total}} = 3T_3 + 2T_2 = \int_0^\tau r_{\text{total}}(X_t)dt$, where $r_{\text{total}} = r_1 + r_2 = (3,2)$

<caption>**Fig. 1.** Phase-type description of the joint distribution of singleton and doubleton branch length ($B_1, B_2$) and total branch length $B_{\text{total}}$. The rewards to the two states correspond to the number of branches of the desired type.</caption>

Intuitively, if the statistic of interest is the number of singleton, doubleton and triplet mutations in a sample of size $n=4$, then a realization of the reward process described here gives the total sums of the branch lengths subtending each of 1,2 or 3 samples in the corresponding genealogical tree. We have previously shown (Hobolth et al., 2019) that the joint branch lengths $B = (B_1, ..., B_d)$ in a time-homogeneous coalescent model are multivariate phase-type distributed. In the infinite sites mutation model, the mutations are sprinkled independently on the various branches according to a Poisson process with mutation rate $\lambda = \theta/2$. Theorem 2.7 in Hobolth et al. (2021) provides the probability generating function for the number of mutations $M = (M_1, ..., M_d)$ in each of the $d$ types of branches with total branch lengths given by $B = (B_1, ..., B_d)$. Specifically, let $B = (B_1, ..., B_d)$ be multivariate phase-type distributed with initial probability vector $\alpha$, sub-intensity matrix $S$, and reward matrix $R$ with non-negative entries. We write $B \sim \text{MPH}(\alpha, S, R)$. If the coalescent process has $p$ states, then the initial state vector is of length $p$, the sub-intensity matrix $S$ is of dimension $p \times p$, and the reward matrix $R$ is of dimension $p \times d$. Furthermore, given $B$, let $M_i$ be independent and Poisson-distributed with rate $\lambda B_i$, i.e. $M_i|B \sim \text{Pois}(\lambda B_i)$, for $i=1, ..., d$.

### Illustrative example: Sample size $n=3$

Consider the number of singletons and doubletons ($M_1, M_2$) in a sample of size $n=3$ as an illustrative example of the phase-type framework for the coalescent process (see Fig. 1). The time $\tau$ to the most recent common ancestor (MRCA) is phase-type distributed with two states. The initial state consists of three singleton branches and the second state consists of one singleton branch and one doubleton branch. Let $\{X_t\}_{0\le t\le \tau}$ denote the state of the coalescent process, i.e. $X_t = 1$ in the times $t$ where the coalescent process has three branches and $X_t = 2$ in the times $t$ where the coalescent process has two branches, and we only consider the ancestral process from the present day $t=0$ to the time $t=\tau$ of the MRCA. The coalescent rate from the first state to the second state is three and the coalescent rate from the second state to the MRCA is one. We write $\tau \sim \text{PH}(\alpha, S)$, where the initial state $\alpha = (1,0)$ and the sub-intensity matrix is given by

$$S = \begin{pmatrix} -3 & 3 \\ 0 & -1 \end{pmatrix}.$$

The total singleton branch length $B_1$ is given by

$$B_1 = \int_0^\tau r_1(X_t)dt,$$

where the reward function $r_1(\cdot)$ is defined by $r_1(1) = 3$ and $r_1(2) = 1$ because state 1 has 3 singleton branches and state 2 has one singleton branch. We summarize the singleton reward function in the vector $r_1 = (3, 1)$. The total doubleton branch length $B_2$ is given by

$$B_2 = \int_0^\tau r_2(X_t)dt,$$

where the reward function $r_2(\cdot)$ is defined by the reward vector $r_2 = (r_2(1), r_2(2)) = (0, 1)$ because state 1 has 0 doubleton branches and state 2 has one doubleton branch. We collect the two reward vectors in a reward matrix $R$ with the two reward vectors as columns.

The marginal distributions of $B_1$ and $B_2$ are so-called reward-transformed phase-type distributions (Section 3.1.8 in Bladt and Nielsen (2017)), and they are both phase-type distributed because the class of phase-type distributions is closed under reward-transformations as long as the rewards are non-negative. The joint distribution of different reward-transformations of the same underlying Markov jump process $\{X_t\}$ is defined to be a multivariate phase-type distribution. The reward-transformations are collected in the matrix $R$ with the reward vectors in the columns. Therefore, the joint distribution $B = (B_1, B_2)$ of the total singleton branch length $B_1$ and the total doubleton branch length $B_2$ is multivariate phase-type distributed $B = (B_1, B_2) \sim \text{MPH}(\alpha, S, R)$ with parameters (initial state, subintensity matrix and reward matrix) given by

$$\alpha = (1,0), \quad S = \begin{pmatrix} -3 & 3 \\ 0 & -1 \end{pmatrix}, \quad \text{and} \quad R = \begin{pmatrix} 3 & 0 \\ 1 & 1 \end{pmatrix}.$$

(1)
Given $\mathbf{B} = (B_1, B_2)$, the number of singleton mutations $M_1$ and number of doubleton mutations $M_2$ are independent Poisson random variables with rates $\lambda B_1$ and $\lambda B_2$, respectively, where $\lambda$ is the mutation rate.♦

In general, the PGF for a stochastic variable $\mathbf{M} = (M_1, \dots, M_d)$ is defined by

$$G(z) = G(z_1, \dots, z_d) = \mathbb{E}[z_1^{M_1} \dots z_d^{M_d}], \quad \text{where } z \in \mathbb{R}^d. \qquad (2)$$

The probability of a genetic sample configuration is determined by differentiating the PGF an appropriate number of times and evaluating at zero argument. The number of times and with which coordinates we need to differentiate is determined by the sample. For non-negative integers $m_1, \dots, m_d$ we have that

$$\mathbb{P}(m_1, \dots, m_d) = (m_1! \dots m_d!)^{-1} \partial_1^{m_1} \dots \partial_d^{m_d} G(\mathbf{0}),$$

where we use the short notation $\mathbb{P}(m_1, \dots, m_d)$ for $\mathbb{P}(M_1 = m_1, \dots, M_d = m_d) = \mathbb{P}(\mathbf{M} = \mathbf{m})$ and $\partial_i^{m_i}$ means derivative with respect to variable $i$ a number $m_i$ times. For example, the probability of a sample configuration $\mathbf{m} = (1, 2, 0, \dots, 0)$ is found by differentiating the PGF once with respect to the first coordinate, twice with respect to the second, evaluating the result at zero, and finally dividing by $1! \cdot 2! = 2$. We refer to Chapter 5 in Grimmett and Stirzaker (1992) for a general introduction to generating functions and examples of their applications in branching processes.

Theorem 2.7 in Hobolth et al. (2021) gives the probability generating function (PGF) for $\mathbf{M} = (M_1, \dots, M_d)$ in the case where $\mathbf{M}$ arises from sprinkling independent Poisson-distributed mutations with rates $\lambda B_1, i = 1, \dots, d$, on multivariate phase-type distributed branches $\mathbf{B} = (B_1, \dots, B_d)$ with initial probability vector $\alpha$, sub-intensity matrix $\mathbf{S}$, and reward matrix $\mathbf{R}$. Let $\mathbf{m} = (m_1, \dots, m_d)$ and define the exit-rate vector $\mathbf{s} = -\mathbf{S}\mathbf{e}$ where $\mathbf{e}$ is the vector with one in every entry. In our illustrative example in Fig. 1 the exit rate vector is $\mathbf{s} = -\mathbf{S}\mathbf{e} = (0, 1)$. Furthermore, let $\Delta(\mathbf{a})$ be the matrix with $\mathbf{a}$ on the diagonal. Then the PGF for $\mathbf{M}$ is given by

$$G(z) = \mathbb{E}[z_1^{M_1} \dots z_d^{M_d}] = \alpha [-\mathbf{S} + \lambda \Delta(\mathbf{R}\mathbf{e}) - \lambda \Delta(\mathbf{R}z)]^{-1} \mathbf{s}. \qquad (3)$$

We define

$$\mathbf{r}_{\text{total}} = \mathbf{R}\mathbf{e}, \qquad (4)$$

to be the total reward vector where each entry is the sum of rewards assigned to the state. In applications $\mathbf{r}_{\text{total}}$ is most often the number of branches in each state. In our illustrative example in Fig. 1 the total number of branches in each of the two states is $\mathbf{r}_{\text{total}} = \mathbf{R}\mathbf{e} = (3, 2)$. We assume the entries in the total reward vector are strictly positive (otherwise we can always remove the states with zero total reward), and re-write the PGF as

$$G(z) = \alpha \underbrace{[I - \lambda^{-1} \Delta(\mathbf{r}_{\text{total}})^{-1} \mathbf{S} - \Delta(\mathbf{r}_{\text{total}})^{-1} \Delta(\mathbf{R}z)]^{-1}}_{\mathbf{A}(z)} \lambda^{-1} \Delta(\mathbf{r}_{\text{total}})^{-1} \mathbf{s} = \alpha \mathbf{A}(z) \lambda^{-1} \Delta(\mathbf{r}_{\text{total}})^{-1} \mathbf{s}. \qquad (5)$$

We emphasize that the generating function is linear in $\mathbf{m}$ before the inversion of the matrix:

$$\mathbf{A}(z) = \underbrace{(I - \lambda^{-1} \Delta(\mathbf{r}_{\text{total}})^{-1} \mathbf{S} - \Delta(\mathbf{r}_{\text{total}})^{-1} \Delta(\mathbf{R}z))}_{\mathbf{C}})^{-1} = (\mathbf{C} - \mathbf{D}(z))^{-1}. \qquad (6)$$

Here, the constant matrix $\mathbf{C}$ is given by

$$\mathbf{C} = I - \lambda^{-1} \Delta(\mathbf{r}_{\text{total}})^{-1} \mathbf{S} = \mathbf{A}(\mathbf{0})^{-1}, \qquad (7)$$

and the diagonal matrix $\mathbf{D}(z)$ is given by

$$\mathbf{D}(z) = \Delta(\mathbf{r}_{\text{total}})^{-1} \Delta(\mathbf{R}z) \qquad (8)$$

with entries that are linear in $z$. This is good news in terms of being able to differentiate with respect to $z$.

In the next few sections we differentiate in general with respect to $\mathbf{m}$ and evaluate at zero. At this point it is instructive to merely consider the probability of zero mutations and describe how it relates to the transition probability matrix for the total number of mutations.

For a sample with zero mutations, the sampling probability equals the PGF evaluated at zero

$$\mathbb{P}(\mathbf{M} = \mathbf{0}) = G(\mathbf{0}) = \alpha \mathbf{A}(\mathbf{0}) \lambda^{-1} \Delta(\mathbf{r}_{\text{total}})^{-1} \mathbf{s}. \qquad (9)$$

From previous work (Theorem 3.5 in Hobolth et al. (2019)) we know that the total number of mutations $M_{\text{total}} = M_1 + \dots + M_d$ is distributed according to the time to absorption of a discrete Markov chain with an absorbing state (or, in other words, is discrete phase-type distributed) with initial distribution $\alpha$ and (sub)probability matrix

$$\mathbf{P}_{\text{total}} = (I - \lambda^{-1} \Delta(\mathbf{r}_{\text{total}})^{-1} \mathbf{S})^{-1} = \mathbf{C}^{-1} = \mathbf{A}(\mathbf{0}). \qquad (10)$$

More precisely, the distribution of the total number of mutations is given by

$$\mathbb{P}(M_{\text{total}} = k) = \alpha \mathbf{P}_{\text{total}}^k p_{\text{total}}, \quad k = 0, 1, \dots, \qquad (11)$$

where

$$p_{\text{total}} = (I - \mathbf{P}_{\text{total}})e \qquad (12)$$

is the exit probability vector.

We have in general that

$$\mathbb{P}(\mathbf{M} = \mathbf{0}) = G(\mathbf{0}) = \alpha \mathbf{A}(\mathbf{0}) \lambda^{-1} \Delta(\mathbf{r}_{\text{total}})^{-1} \mathbf{s} = \alpha p_{\text{total}}. \qquad (13)$$
This result follows from (9) and the fact that $\mathbf{A}(0)\lambda^{-1} \triangle (\mathbf{r}_{\text{total}})^{-1}s = p_{\text{total}}$. To prove the latter result we first note that from the general matrix equality $(\mathbf{A} - \mathbf{B})^{-1}(\mathbf{A} - \mathbf{B}) = \mathbf{I}$ we get $(\mathbf{A} - \mathbf{B})^{-1}\mathbf{A} - (\mathbf{A} - \mathbf{B})^{-1}\mathbf{B} = \mathbf{I}$. Second we make a direct computation to show the desired result:

$$ 
\begin{aligned}
\mathbf{A}(0)\lambda^{-1} \triangle (\mathbf{r}_{\text{total}})^{-1}s &= \left(\mathbf{I} - \lambda^{-1} \triangle (\mathbf{r}_{\text{total}})^{-1}\mathbf{S}\right)^{-1} \lambda^{-1} \triangle (\mathbf{r}_{\text{total}})^{-1}(-\mathbf{S}e) \\
&= -\left(\mathbf{I} - \lambda^{-1} \triangle (\mathbf{r}_{\text{total}})^{-1}\mathbf{S}\right)^{-1} \lambda^{-1} \triangle (\mathbf{r}_{\text{total}})^{-1}\mathbf{S}e \\
&= \left[\mathbf{I} - \left(\mathbf{I} - \lambda^{-1} \triangle (\mathbf{r}_{\text{total}})^{-1}\mathbf{S}\right)^{-1}\right]e \\
&= (\mathbf{I} - \mathbf{P}_{\text{total}})e \\
&= p_{\text{total}}.
\end{aligned}
\quad (14) $$

### Illustrative example: Sample size $n = 3$ (continued)

In our illustrative case (1) of the number of singletons and doubletons we find from (7) that

$$ C = \mathbf{I} - \lambda^{-1} \triangle (\mathbf{r}_{\text{total}})^{-1}\mathbf{S} = \begin{pmatrix} 1 + \frac{1}{\lambda} & -\frac{1}{\lambda} \\ 0 & 1 + \frac{1}{2\lambda} \end{pmatrix} = \begin{pmatrix} \frac{(1 + \lambda)}{\lambda} & -\frac{1}{\lambda} \\ 0 & \frac{(1 + 2\lambda)}{2\lambda} \end{pmatrix}, $$

which means that

$$ \mathbf{P}_{\text{total}} = \mathbf{C}^{-1} = \begin{pmatrix} \frac{3\lambda}{(3 + 3\lambda)} & \frac{3}{(3 + 3\lambda)}\frac{2\lambda}{(1 + 2\lambda)} \\ 0 & \frac{2\lambda}{(1 + 2\lambda)} \end{pmatrix}. \quad (15) $$

We emphasize the nice interpretation of $\mathbf{P}_{\text{total}}$ as a transition probability matrix for a mutation at a level in the tree when three lineages are present (state 1) or when two lineages are present (state 2): a transition $1 \to 1$ corresponds to a mutation (rate $3\lambda$) before a coalescence (rate 3) when three lineages are present; a transition $1 \to 2$ corresponds to a coalescence when three lineages are present followed by a mutation (rate $2\lambda$) before a coalescence (rate 1) when two lineages are present; a transition $2 \to 1$ is not possible; a transition $2 \to 2$ corresponds to a mutation before a coalescence when two lineages are present.

Furthermore, we find the exit probability vector

$$ \mathbf{p}_{\text{total}} = (\mathbf{I} - \mathbf{P}_{\text{total}})e = \begin{pmatrix} \frac{3}{(3+3\lambda)} & \frac{1}{(1+2\lambda)} \\ \frac{1}{(1+2\lambda)} & \frac{1}{(1+2\lambda)} \end{pmatrix} \quad (16) $$

which is the probability of having no more mutations starting from either three lineages (first entry) or two lineages (second entry). The interpretation of the expression (11) for the total number of mutations is now clear: $\alpha$ is the probability vector for beginning in the various states, $\mathbf{P}_{\text{total}}$ is the transition probability matrix for a mutation, and $\mathbf{p}_{\text{total}}$ is the probability vector for exiting the system from the different states.

We find from (13) that the probability of zero mutations is

$$ G(0, 0) = \alpha \mathbf{p}_{\text{total}} = \frac{3}{(3 + 3\lambda)} \frac{1}{(1 + 2\lambda)}, $$

and this expression corresponds to the two coalescent events taking place before any mutation events. ◦

In Section 4 we differentiate the PGF in general with respect to the variables $(m_1, \dots, m_d)$ an appropriate number of times with respect to each variable and evaluate at zero. Fortunately, for time-homogeneous coalescent models, this task can be accomplished using a matrix-analytical formula, as we will demonstrate in Section 4.1. However, before introducing the general sampling formula we provide in the next section a number of illustrative examples for a detailed understanding of the PGF.

## 3. Four illustrative examples for understanding the PGF

The benefit of the PGF for deriving concise expressions regarding reward processes may seem opaque, as the mathematical manipulations are done in the transformed space (and thus obscured by the ‘Laplacian Curtain’, see e.g. the discussion in the preface of Neuts (1981)). We now give four simple examples that we find useful for getting familiar with the notation and for obtaining intuition about the PGF. In Section 5 we give more examples, and in Section 6 we provide a more involved application in order to illustrate the easy implementation and powerful machinery of the matrix-analytical framework.

### 3.1. Example 1: Total number of mutations for sample size $n = 2$

Consider the number of mutations $M$ for a sample of size two. We have that the time to the most recent common ancestor $B$ is exponentially distributed with rate one, so $\alpha = 1$, $S = -1$, and $s = 1$. The number of lineages is two, so we give a reward of $R = 2$ to the state. Inserting in (3) we get the PGF

$$ G(z) = 1 \cdot (1 + 2\lambda - 2\lambda z)^{-1} \cdot 1 = \frac{1}{1 + 2\lambda - 2\lambda z}, \quad \text{for } z \in \mathbb{R}. $$

It is now easy to see that

$$ \mathbb{P}(M = 0) = G(0) = \frac{1}{1 + 2\lambda}, \quad \mathbb{P}(M = 1) = \partial^1 G(0) = \frac{2\lambda}{(1 + 2\lambda)^2} = \frac{2\lambda}{(1 + 2\lambda)} \frac{1}{(1 + 2\lambda)}, $$

and in general for $m \ge 1$ we have

$$ \mathbb{P}(M = m) = \frac{1}{m!} \partial^m G(0) = \frac{m!}{m! (1 + 2\lambda)^{m+1}} = \left(\frac{2\lambda}{1 + 2\lambda}\right)^m \frac{1}{(1 + 2\lambda)}. $$
We see that the number of mutations is geometrically distributed with probability of success (i.e. coalescent) $1/(1+2\lambda)$ and probability of failure (i.e. mutation) $2\lambda/(1+2\lambda)$. This is of course a consequence of the two competing events and their rates: coalescent with a rate of one and mutation with a rate of $2\lambda$.

### 3.2. Example 2: Number of mutations in each branch for sample size $n = 2$

Consider the number of mutations ($M_1, M_2$) in each branch for a sample of size two. The time to the most recent common ancestor $B$ is again exponentially distributed with rate one, so $\alpha = 1$, $S = -1$, and $s = 1$. We distinguish between the two branches, and therefore the reward vector is $\mathbf{R} = (1, 1)$. From (4) we get $\mathbf{r}_{\text{total}} = \Delta(\mathbf{R}\mathbf{e}) = \mathbf{R}\mathbf{e} = 2$, from (7) we get

$$C = 1 + \frac{1}{2\lambda} = \frac{1 + 2\lambda}{2\lambda},$$

and from (8) we get

$$D(z_1, z_2) = \frac{z_1 + z_2}{2}.$$

Inserting in firstly (6) and secondly (5) we get the PGF

$$G(z_1, z_2) = 1 \cdot \left( \frac{1 + 2\lambda}{2\lambda} - \frac{z_1 + z_2}{2} \right)^{-1} \cdot \frac{1}{2\lambda} \cdot 1 = \left( \frac{1 + 2\lambda}{2\lambda} - \frac{z_1 + z_2}{2} \right)^{-1} \frac{1}{2\lambda}, \quad \text{for } (z_1, z_2) \in \mathbb{R}^2.$$

We see that

$$P(M_1 = 0, M_2 = 0) = G(0, 0) = \frac{2\lambda}{(1 + 2\lambda) 2\lambda} = \frac{1}{1 + 2\lambda}.$$

Furthermore we have

$$\begin{aligned}
P(M_1 = 1, M_2 = 0) &= \partial_1^1 G(0, 0) = \left( \frac{1 + 2\lambda}{2\lambda} - \frac{z_1 + z_2}{2} \right)^{-2} \frac{1}{2\lambda} \Big|_{(z_1, z_2)=(0,0)} \\
&= \frac{2\lambda}{(1 + 2\lambda)^2} \frac{1}{2} = \frac{2\lambda}{(1 + 2\lambda) \frac{1}{2}} = \frac{1}{(1 + 2\lambda)},
\end{aligned}$$

which corresponds to a mutation before the coalescent (probability $2\lambda/(1+2\lambda)$), the mutation must occur in the first branch (probability $1/2$), and finally a coalescent (probability $1/(1+2\lambda)$).

In general for non-negative integers $(m_1, m_2)$ we get

$$\begin{aligned}
P(M_1 = m_1, M_2 = m_2) &= (m_1! m_2!)^{-1} \partial_1^{m_1} \partial_2^{m_2} G(0, 0) \\
&= \frac{(m_1 + m_2)!}{m_1! m_2!} \left\{ \frac{2\lambda}{1 + 2\lambda} \right\}^{m_1 + m_2 + 1} \left(\frac{1}{2}\right)^{m_1} \left(\frac{1}{2}\right)^{m_2} \frac{1}{2\lambda} \\
&= \binom{m_1 + m_2}{m_1} \left(\frac{1}{2}\right)^{m_1 + m_2} \left\{ \frac{2\lambda}{1 + 2\lambda} \right\}^{m_1 + m_2} \frac{1}{(1 + 2\lambda)}.
\end{aligned}$$

This expression corresponds to a binomial distribution with $m_1 + m_2$ trials (number of mutations) with equal probability for the two cases, and a geometric distribution for the total number of mutations with probability of mutation $2\lambda/(1+2\lambda)$.

Actually, we can observe that the total number of mutations is discrete phase-type distributed with initial distribution $\alpha = 1$, (sub)probability (10) given by $P_{\text{total}} = C^{-1} = (2\lambda)/(1+2\lambda)$ and exit probability vector (12) given by $p_{\text{total}} = (1 - P_{\text{total}}) = 1/(1+2\lambda)$. These observations are in accordance with the previous example because a discrete phase-type distribution of dimension one is a geometric distribution or, equivalently, the time to absorption in a discrete Markov chain of dimension two where one of the states is an absorbing state.

### 3.3. Example 3: Site frequency spectrum for sample size $n = 3$

This is our illustrative example from Fig. 1 in Section 2, and in particular the initial state vector $\alpha$, the sub-intensity matrix $S$ and the reward matrix $\mathbf{R}$ are given by (1). From (8) we get

$$\mathbf{D}(z_1, z_2) = \Delta(\mathbf{r}_{\text{total}})^{-1} \Delta(\mathbf{R}\mathbf{z}) = \begin{pmatrix} 1/3 & 0 \\ 0 & 1/2 \end{pmatrix} \begin{pmatrix} 3z_1 & 0 \\ 0 & z_1 + z_2 \end{pmatrix} = \begin{pmatrix} z_1 & 0 \\ 0 & (z_1 + z_2)/2 \end{pmatrix}.$$

In order to determine the probability of one singleton and zero doubleton mutations we need to calculate the first derivative of the PGF with respect to $m_1$, and evaluate at $(m_1, m_2) = (0, 0)$. Recall that the derivative of an inverse matrix is

$$\partial(\mathbf{U}(\mathbf{z})^{-1}) = -\mathbf{U}(\mathbf{z})^{-1} (\partial\mathbf{U}(\mathbf{z})) \mathbf{U}(\mathbf{z})^{-1},$$

and we find (since $\mathbf{A}(\mathbf{z})^{-1} = \mathbf{C} - \mathbf{D}(\mathbf{z})$) that

$$\partial_1 \mathbf{A}(\mathbf{z}) = \mathbf{A}(\mathbf{z}) (\partial_1 \mathbf{D}(\mathbf{z})) \mathbf{A}(\mathbf{z}),$$

where the notation $\partial_1$ means derivative with respect to $m_1$. We see from (8) that

$$\partial_1 \mathbf{D}(\mathbf{z}) = \Delta(\mathbf{r}_{\text{total}})^{-1} \Delta(\mathbf{R}\mathbf{e}_1) = \Delta(\mathbf{r}_{\text{total}})^{-1} \Delta(\mathbf{r}_1) = \Delta(\mathbf{r}_1/\mathbf{r}_{\text{total}}),$$

where vector division is element-wise, $e_1 = (1, 0)$, and $\mathbf{r}_1 = (3, 1)$ is the vector of rewards for singletons. We define $\mathbf{D}_1 = \Delta(\mathbf{r}_1/\mathbf{r}_{\text{total}})$ and note that the vector $\mathbf{r}_1/\mathbf{r}_{\text{total}} = (1, 1/2)$ is the vector of fractions of singleton branches out of the total number of branches in each of the two states. From (5)
and (14) we get

$$ 
\begin{aligned}
\mathbb{P}(M_1 = 1, M_2 = 0) &= \partial_1 G(0, 0) = \alpha \mathbf{A}(0) \mathbf{D}_1 \mathbf{A}(0) \lambda^{-1} \Delta (\mathbf{r}_{\text{total}})^{-1} s = \alpha \mathbf{P}_{\text{total}} \mathbf{D}_1 \mathbf{P}_{\text{total}} \\
&= \frac{3\lambda}{(3+3\lambda)} \frac{3}{(3+3\lambda)} \frac{1}{(1+2\lambda)} + \frac{3}{(3+3\lambda)} \frac{\lambda}{(1+2\lambda)} \frac{1}{(1+2\lambda)},
\end{aligned}
$$

and this expression corresponds to either a mutation in the ancestral process during the time with three lineages (the first term) or a mutation in the singleton branch during the time with two lineages (the second term).

We can calculate the probability of zero singleton mutations and one doubleton mutation using a similar procedure as above. The doubleton reward vector is $\mathbf{r}_2 = (0, 1)$ and we get $\mathbf{r}_2/\mathbf{r}_{\text{total}} = (0, 1/2)$; in the first state all three branches are singletons, and the second state consists of one singleton branch and one doubleton branch. We define $\mathbf{D}_2 = \Delta (\mathbf{r}_2/\mathbf{r}_{\text{total}})$ and get

$$ 
\begin{aligned}
\mathbb{P}(M_1 = 0, M_2 = 1) &= \partial_2 G(0, 0) = \alpha \mathbf{A}(0) \mathbf{D}_2 \mathbf{A}(0) \lambda^{-1} \Delta (\mathbf{r}_{\text{total}})^{-1} s = \alpha \mathbf{P}_{\text{total}} \mathbf{D}_2 \mathbf{P}_{\text{total}} \\
&= \frac{3}{(3+3\lambda)} \frac{\lambda}{(1+2\lambda)} \frac{1}{(1+2\lambda)},
\end{aligned}
$$

and this expression corresponds to a single mutation in the doubleton branch.

In order to calculate the probability for a sample with 2 singleton mutations and 1 doubleton mutation we would need to calculate

$$ \frac{1}{2} \partial_1^2 \partial_2^1 \mathbf{A}(0) = (\mathbf{A}(0) \mathbf{D}_1)^2 (\mathbf{A}(0) \mathbf{D}_2) + (\mathbf{A}(0) \mathbf{D}_1)(\mathbf{A}(0) \mathbf{D}_2)(\mathbf{A}(0) \mathbf{D}_1) + (\mathbf{A}(0) \mathbf{D}_2)(\mathbf{A}(0) \mathbf{D}_1)^2, \quad (17) $$

and multiply by $\alpha$ from the left and $\mathbf{A}(0)\lambda^{-1} \Delta (\mathbf{r}_{\text{total}})^{-1} s = \mathbf{p}_{\text{total}}$ from the right. If we let $\mathbf{P}_1 = \mathbf{A}(0)\mathbf{D}_1 = \mathbf{P}_{\text{total}}\mathbf{D}_1$ and $\mathbf{P}_2 = \mathbf{A}(0)\mathbf{D}_2 = \mathbf{P}_{\text{total}}\mathbf{D}_2$ we have

$$ \mathbb{P}(M_1 = 2, M_2 = 1) = \alpha (\mathbf{P}_1 \mathbf{P}_1 \mathbf{P}_2 + \mathbf{P}_1 \mathbf{P}_2 \mathbf{P}_1 + \mathbf{P}_2 \mathbf{P}_1 \mathbf{P}_1) \mathbf{p}_{\text{total}}. \quad (18) $$

Similarly, in order to calculate the probability of 2 singleton mutations and 2 doubleton mutations we would need to calculate the mixed higher-order derivative

$$ \frac{1}{4} \partial_1^2 \partial_2^2 \mathbf{A}(0) = \mathbf{P}_1^2 \mathbf{P}_2^2 + \mathbf{P}_1 \mathbf{P}_2 \mathbf{P}_1 \mathbf{P}_2 + \mathbf{P}_1 \mathbf{P}_2^2 \mathbf{P}_1 + \mathbf{P}_2 \mathbf{P}_1^2 \mathbf{P}_2 + \mathbf{P}_2 \mathbf{P}_1 \mathbf{P}_2 \mathbf{P}_1 + \mathbf{P}_2^2 \mathbf{P}_1^2, \quad (19) $$

and multiply by $\alpha$ from the left and $\mathbf{p}_{\text{total}}$ from the right:

$$ \mathbb{P}(M_1 = 2, M_2 = 2) = \alpha (\mathbf{P}_1^2 \mathbf{P}_2^2 + \mathbf{P}_1 \mathbf{P}_2 \mathbf{P}_1 \mathbf{P}_2 + \mathbf{P}_1 \mathbf{P}_2^2 \mathbf{P}_1 + \mathbf{P}_2 \mathbf{P}_1^2 \mathbf{P}_2 + \mathbf{P}_2 \mathbf{P}_1 \mathbf{P}_2 \mathbf{P}_1 + \mathbf{P}_2^2 \mathbf{P}_1^2) \mathbf{p}_{\text{total}}. \quad (20) $$

The sampling formulas involve products of the matrices $\mathbf{P}_i = \mathbf{P}_{\text{total}}\mathbf{D}_i = \mathbf{P}_{\text{total}} \Delta (\mathbf{r}_i/\mathbf{r}_{\text{total}})$, where the fraction $\mathbf{r}_i/\mathbf{r}_{\text{total}}$ of the two vectors should be calculated element-wise. It is natural to view the expressions in terms of a thinning process of the total number of mutations, where the thinning is state-dependent, and for each state j the probability that a mutation is of type i is given by $\mathbf{r}_i(j)/\mathbf{r}_{\text{total}}(j)$. In this example we find that

$$ \mathbf{P}_1 = \mathbf{P}_{\text{total}} \mathbf{D}_1 = \mathbf{P}_{\text{total}} \Delta (\mathbf{r}_1 / \mathbf{r}_{\text{total}}) = \begin{pmatrix} \frac{3\lambda}{3+3\lambda} & \frac{3}{(3+3\lambda)} & \frac{\lambda}{(1+2\lambda)} \\ 0 & \frac{\lambda}{(1+2\lambda)} & 0 \end{pmatrix} \quad (21) $$

and

$$ \mathbf{P}_2 = \mathbf{P}_{\text{total}} \mathbf{D}_2 = \mathbf{P}_{\text{total}} \Delta (\mathbf{r}_2 / \mathbf{r}_{\text{total}}) = \begin{pmatrix} 0 & \frac{3}{(3+3\lambda)} & \frac{\lambda}{(1+2\lambda)} \\ 0 & \frac{\lambda}{(1+2\lambda)} & 0 \end{pmatrix}. \quad (22) $$

We note that because the matrices $\mathbf{P}_i$ correspond to a mutation of type i we have that a term such as $\alpha \mathbf{P}_1 \mathbf{P}_1 \mathbf{P}_2 \mathbf{P}_{\text{total}}$ equals the probability of exactly three mutations in the ancestral process where the two mutations of type one (the two singleton mutations) occur before the single mutation of type two (the doubleton mutation). In (18) we find the probability of two singleton mutations and one doubleton mutation by summing over all three ordered ways that these three mutations can arise.

We find it useful and easy to calculate derivatives of the type (19) by using the `multiset()` command in the `partitions` package in R (Hankin and West, 2007). For example the command

`> multiset(c(1, 1, 2, 2))`

returns the matrix

[1,] 1 1 1 2 2 2
[2,] 1 2 2 1 1 2
[3,] 2 1 2 1 2 1
[4,] 2 2 1 2 1 1

The columns in this matrix corresponds to all possible ordered ways that two mutations of type 1 and two mutations of type 2 can arise. This makes it easy to implement a general function for calculating the desired derivatives: a loop across rows where the corresponding $\mathbf{P}_i$'s are multiplied and a loop across columns where the products are added corresponds to the expression in (19). The fact that the matrices $\mathbf{P}_i = \mathbf{P}_{\text{total}}\mathbf{D}_i$, $i = 1, 2$, are sub-transition probability matrices makes the evaluation computationally very stable because all entries in the product of the matrices are between zero and one.

We can summarize the procedure for calculating the probability of a site frequency spectrum ($M_1, M_2$) for three sequences in the following four simple steps.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Phase-type description of the joint distribution of singleton, doubleton and tripleton branch lengths ($B_1, B_2, B_3$). The rewards to each of the four states correspond to the number of branches of the desired type.</figcaption></figure>

1. Define the evolutionary model in terms of the time-homogeneous coalescent process ($\alpha, S, R$) from Eq. (1) and the mutation rate $\lambda$.

2. Calculate $r_{\text{total}} = Re, P_{\text{total}}$ from (10), $p_{\text{total}}$ from (16), $D_i = \Delta(r_i/r_{\text{total}})$, $i = 1, 2$, and $P_i = P_{\text{total}} D_i$, $i = 1, 2$. The resulting (sub)transition matrices are given in (21) and (22).

3. Use the multiset() function in the R package partitions to enumerate all possible ways that a sample ($m_1, m_2$) of $m_1$ singletons and $m_2$ doubletons may arise, and calculate the corresponding matrix products and additions. For example, if the sample is ($m_1, m_2$) = (2, 2), then calculate (19).

4. Multiply the resulting matrix by $\alpha$ from the left and $p_{\text{total}}$ from the right; the final result is the probability of the sample. For example, if the sample is ($m_1, m_2$) = (2, 2), then calculate (20).

In the next section we generalize this procedure to calculate the probability of a sample from a general time-homogeneous coalescent model.

### 3.4. Example 4: site frequency spectrum for sample size $n = 4$

In this case we have four states: a state with four singleton branches, a state with two singleton branches and one doubleton branch, a state with one singleton and one tripleton branch, and a state with two doubleton branches. See Fig. 2 for an overview and enumeration of the states.

The joint distribution of the total singleton branch length $B_1$, total doubleton branch length $B_2$, and total tripleton branch length $B_3$ is multivariate phase-type distributed $B = (B_1, B_2, B_3) \sim \text{MPH}(\alpha, S, R)$ with beginning state $\alpha$, (sub)intensity matrix $S$ and reward matrix $R$ given by

$$ \alpha = (1, 0, 0, 0), \quad S = \begin{pmatrix} -6 & 6 & 0 & 0 \\ 0 & -3 & 2 & 1 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix} \quad \text{and} \quad R = \begin{pmatrix} 4 & 0 & 0 \\ 2 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 2 & 0 \end{pmatrix}. $$

We find that

$$ r_{\text{total}} = Re = (4, 3, 2, 2), $$

which corresponds to the number of branches in each state. Using Eq. (10) we get

$$ \begin{aligned} P_{\text{total}} &= \begin{pmatrix} 1 + \frac{1}{\lambda} & \frac{1}{\lambda} & 0 & 0 \\ 0 & 1 + \frac{1}{\lambda} & -\frac{1}{\lambda} & \frac{1}{3} \\ 0 & 0 & 1 + \frac{1}{\lambda} & 0 \\ 0 & 0 & 0 & 1 + \frac{1}{\lambda} \end{pmatrix}^{-1} \\ &= \begin{pmatrix} \frac{4\lambda}{(6+4\lambda)} & \frac{6}{(6+4\lambda)} & \frac{3\lambda}{(3+3\lambda)} & \frac{6}{(6+4\lambda)} & \frac{3}{(3+3\lambda)} & \frac{2}{3} & \frac{2\lambda}{(1+2\lambda)} & \frac{6}{(6+4\lambda)} & \frac{3}{(3+3\lambda)} & \frac{1}{3} & \frac{2\lambda}{(1+2\lambda)} \\ 0 & \frac{3\lambda}{(3+3\lambda)} & \frac{2}{3} & \frac{2\lambda}{(1+2\lambda)} & \frac{3}{(3+3\lambda)} & \frac{1}{3} & \frac{2\lambda}{(1+2\lambda)} & \frac{3}{(3+3\lambda)} & \frac{1}{3} & \frac{2\lambda}{(1+2\lambda)} & \\ 0 & 0 & \frac{2\lambda}{(1+2\lambda)} & 0 & \frac{2\lambda}{(1+2\lambda)} & 0 & 0 & \frac{2\lambda}{(1+2\lambda)} & 0 & 0 \end{pmatrix} \end{aligned} $$

and from Eq. (12) we get

$$ P_{\text{total}} = \left( \frac{6}{(6+4\lambda)} \frac{3}{(3+3\lambda)} \frac{1}{(1+2\lambda)}, \frac{3}{(3+3\lambda)} \frac{1}{(1+2\lambda)}, \frac{1}{(1+2\lambda)}, \frac{1}{(1+2\lambda)} \right). $$

We note that the entries in $p_{\text{total}}$ corresponds to finding common ancestry with no further mutations from state 1, 2, 3 or 4, respectively.
We have the thinning fractions in each of the four states for each of the three mutation types

$$r_1/r_{\text{total}} = (1, 2/3, 1/2, 0), \quad r_2/r_{\text{total}} = (0, 1/3, 0, 1), \quad r_3/r_{\text{total}} = (0, 0, 1/2, 0),$$

and we find the sub-transition probabilities for each of the three mutation types (singleton, doubletons or tripletons) to be

$$P_1 = P_{\text{total}} \Delta (r_1/r_{\text{total}}) = \begin{pmatrix} \frac{4\lambda}{(6+4\lambda)} & \frac{6}{(6+4\lambda)} & \frac{2\lambda}{(3+3\lambda)} & \frac{6}{(6+4\lambda)} & \frac{3}{(3+3\lambda)} & \frac{2}{3} & \frac{\lambda}{(1+2\lambda)} & 0 \\ 0 & \frac{2\lambda}{(3+3\lambda)} & 0 & \frac{3}{(3+3\lambda)} & \frac{2}{3} & \frac{\lambda}{(1+2\lambda)} & 0 & 0 \\ 0 & 0 & 0 & \frac{\lambda}{(1+2\lambda)} & 0 & 0 & 0 & 0 \end{pmatrix},$$

$$P_2 = P_{\text{total}} \Delta (r_2/r_{\text{total}}) = \begin{pmatrix} 0 & \frac{6}{(6+4\lambda)} & \frac{\lambda}{(3+3\lambda)} & 0 & \frac{6}{(6+4\lambda)} & \frac{3}{(3+3\lambda)} & \frac{1}{3} & \frac{2\lambda}{(1+2\lambda)} \\ 0 & \frac{\lambda}{(3+3\lambda)} & 0 & 0 & \frac{3}{(3+3\lambda)} & \frac{1}{3} & \frac{2\lambda}{(1+2\lambda)} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{2\lambda}{(1+2\lambda)} & 0 & 0 & 0 \end{pmatrix},$$

and

$$P_3 = P_{\text{total}} \Delta (r_3/r_{\text{total}}) = \begin{pmatrix} 0 & 0 & \frac{6}{(6+4\lambda)} & \frac{3}{(3+3\lambda)} & \frac{2}{3} & \frac{\lambda}{(1+2\lambda)} & 0 \\ 0 & 0 & \frac{3}{(3+3\lambda)} & \frac{2}{3} & \frac{\lambda}{(1+2\lambda)} & 0 & 0 \\ 0 & 0 & \frac{\lambda}{(1+2\lambda)} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}.$$

We emphasize that the interpretation of the transition (sub)probability matrices is rather clear. For example, we see that $P_3P_2 = 0$ because doubleton mutations are not possible after tripleton mutations, the last column in both $P_1$ and $P_3$ is the zero vector because singleton and tripleton mutations are not possible in state 4, and terms such as $2\lambda/(3+3\lambda)$ indicate a mutation in two out of three branches (e.g. a singleton mutation in state 2).

A general construction of the beginning state $\alpha$, sub-intensity matrix $S$ and reward matrix $R$ for the block-counting process of the Kingman coalescent for a general sample size $n$ is described in Algorithm 4.2 in Hobolth et al. (2019), and the corresponding implementation is available in our GitHub repository.

## 4. General sampling formula and algorithm

In this section, we describe the general sampling formula and corresponding algorithm for calculating the probability of a sample. In order to describe the sampling formula we need to introduce some mathematical notation for all permutations of a list.

Recall that the `multiset()` function from the R package `partitions` (Hankin and West, 2007) allows for an easy enumeration of all possible ways that the mutations may arise in the ancestral process. In mathematical notation, let $\sigma = \sigma(m_1, \dots, m_d)$ be the matrix of all permutations of the multiset of the configuration

$$\underbrace{(1, \dots, 1)}_{m_1 \text{ repeats}} \underbrace{2, \dots, 2}_{m_2 \text{ repeats}} \dots \underbrace{d, \dots, d}_{m_d \text{ repeats}}$$

For example

$$\sigma(2, 1, 1) = \begin{pmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 2 & 2 & 2 & 3 & 3 & 3 \\ 1 & 1 & 2 & 2 & 3 & 3 & 1 & 1 & 3 & 1 & 1 & 2 \\ 2 & 3 & 1 & 3 & 1 & 2 & 1 & 3 & 1 & 1 & 2 & 1 \\ 3 & 2 & 3 & 1 & 2 & 1 & 3 & 1 & 1 & 2 & 1 & 1 \end{pmatrix}.$$

The number of mutations (rows of $\sigma$) is $r_{\sigma} = m_{\text{total}} = m_1 + \dots + m_d$, and the number of permutations (columns of $\sigma$) is the multinomial coefficient $c_{\sigma} = \binom{m_{\text{total}}}{m_1 \dots m_d}$. In this example we have $r_{\sigma} = 2+1+1=4$ and $c_{\sigma} = \binom{4}{2} = 12$. The sampling formula is now given in the following theorem.

### 4.1. Theorem: Sampling formula for the infinite sites mutation model

Consider a time-homogeneous coalescent process with $p$ states and $d$ types of branches for tracking mutations. The coalescent process is defined by an initial probability vector $\alpha$ of length $p$, a $p \times p$ sub-intensity matrix $S$, and a $p \times d$ non-negative reward matrix $R$ with columns $r_i$, $i = 1, \dots, d$. Assume the infinite sites mutation model with a mutation rate $\lambda > 0$. Let $r_{\text{total}} = Re$ and assume all entries in $r_{\text{total}}$ are positive. Let

$$P_{\text{total}} = (I - \lambda^{-1} \Delta (r_{\text{total}})^{-1} S)^{-1}, \quad P_{\text{total}} = (I - P_{\text{total}})e, \\ D_i = \Delta (r_i / r_{\text{total}}) \quad \text{and} \quad P_i = P_{\text{total}} D_i, \quad i = 1, \dots, d.$$
Let $\sigma = \sigma(m)$ be the permutation matrix of the $d$-dimensional sample of non-negative integers $m \in \mathbb{N}_{\ge 0}^d$. The permutation matrix has entries $\sigma_{ij}$, number of rows $r_\sigma$, and number of columns $c_\sigma$. The probability of the sample is given by

$$P(m) = \alpha \left( \sum_{j=1}^{c_\sigma} \prod_{i=1}^{r_\sigma} P_{\sigma_{ij}} \right) p_{\text{total}}$$

**Proof.** The theorem follows from taking an appropriate number of derivatives of the probability generating function (5) and evaluate at zero. This line of reasoning was demonstrated in detail in Example 3 in the previous section. We first take advantage of the linearity of the diagonal matrix $D(z)$ and find

$$\partial_i \{ D(z) \} = \Delta (r_i / r_{\text{total}}) = D_i,$$

and second use calculus for taking partial derivatives of an inverse matrix to get

$$\partial_i \{ A(z) \} = A(z) D_i A(z).$$

Recall from (10) that $A(0) = P_{\text{total}}$, and using (14) we get the probability $\alpha P_{\text{total}} D_i P_{\text{total}} = \alpha P_i P_{\text{total}}$ for a sample with one mutation of type $i$ and zero mutations of all other types.

In order to calculate the probability of one mutation of type $i$, one mutation of type $j$, and zero mutations of all other types we find

$$\partial_j \{ \partial_i \{ A(z) \} \} = A(z) D_j A(z) D_i A(z) + A(z) D_i A(z) D_j A(z),$$

and this translates into the probability $\alpha(P_i P_j + P_j P_i) p_{\text{total}}$.

Applying these formulas recursively allows us to formulate a general expression of $\partial_1^{m_1} \cdots \partial_d^{m_d} \{A(z)\}$ in terms of the multiset notation, and this expression can then be turned into the sampling formula by evaluating in zero and multiplying by the appropriate terms on the left side and right side of the expression.

Basically, the theorem is an application of the Faà de Bruno formula (see e.g. Hardy, 2006). □

The sampling formula for calculating $P(m)$ consists of $c_\sigma$ summands, and for each summand we have $(r_\sigma - 1)$ matrix multiplications. This results in a total number of $c_\sigma (r_\sigma - 1)$ matrix multiplications. We refer to Section 7.1 for a more elaborate discussion about computational complexity. In particular we discuss the number of summands $c_\sigma$ and how the computational complexity can potentially be reduced by pre-computing matrix powers.

## 4.2. Algorithm: Probability of a genetic sample

The sampling formula can be implemented in an algorithm consisting of four simple steps.

### Algorithm: Probability of a genetic sample

1. **Input:** Evolutionary model ($\alpha$, $S$, $R$, $\lambda$) and sample $m$.

The coalescent process consists of three ingredients: an initial probability vector $\alpha$ of length $p$, a $p \times p$ sub-intensity matrix $S$, and a $p \times d$ non-negative reward matrix $R$. Mutations arise according to the infinite sites model with rate $\lambda$. The sample $m$ is a $d$-dimensional vector.

2. **Calculate**

$$r_{\text{total}} = Re, \quad P_{\text{total}} = (I - \lambda^{-1} \Delta (r_{\text{total}})^{-1} S)^{-1}, \quad p_{\text{total}} = (I - P_{\text{total}})e, \\ D_i = \Delta (r_i / r_{\text{total}}), \quad P_i = P_{\text{total}} D_i, \quad i = 1, \dots, d.$$

3. Determine the permutation matrix $\sigma(m)$ of the sample and calculate the matrix $\sum_{j=1}^{c_\sigma} \prod_{i=1}^{r_\sigma} P_{\sigma_{ij}}$, where $r_\sigma$ is the number of rows and $c_\sigma$ is the number of columns of $\sigma(m)$.

4. Multiply the matrix from Step 3 by $\alpha$ from the left and $p_{\text{total}}$ from the right to obtain the probability of the sample, i.e. calculate $P(m) = \alpha (\sum_{j=1}^{c_\sigma} \prod_{i=1}^{r_\sigma} P_{\sigma_{ij}}) p_{\text{total}}$.

**Output:** Probability of the sample $P(m)$.

In the next two sections we provide several examples of the algorithm. It is important to emphasize that the construction of the sub-intensity matrix $S$ and reward matrix $R$ in Step 1 is the only difficult part of the algorithm. Fortunately, most often the construction of $S$ and $R$ can be automatized. For example, to calculate the probability of a site frequency spectrum for any given sample size $n$, the two matrices are given by the block-counting process, and a function to determine the two matrices for this process is available in PhaseTypeR (Rivas-González et al., 2023). A similar general construction is available in Section 6.2 for the Beta-coalescent and can be found in our GitHub repository at https://github.com/AsgerHobolth/MatrixAnalyticalSamplingFormula.

The size of the state space for a sample size $n$ for the standard block counting process equals the partition function $p(n)$ from number theory (e.g. Table 1 in Hobolth et al. (2024)). For sample sizes $n = (2, 4, 6, 10, 15, 20, 25, 30)$ the partition function is equal to $p(n) = (2, 5, 11, 42, 176, 627, 1958, 5,604)$. Our framework can thus be used to calculate exact probabilities for site frequency spectra with sample sizes $n \le 25$.

Exact calculation of the joint allele frequency spectrum is investigated for the structured coalescent model with 2 islands in Kern and Hey (2017). If the sample sizes $(n_1, n_2)$ in the two islands are equal $(n_1 = n_2 = n/2)$ they get for total sample size $n = (2, 4, 6)$ that the size of the state space equals $(6, 46, 268)$. For $n = 12$ $(n_1 = n_2 = 6)$ the state space can be seen from Figure 2 in Kern and Hey (2017) to be larger than $\exp(10)$, which is larger than 22.000. Kern and Hey (2017) mention that their framework can handle sample sizes $n_1 = n_2 < 9$. A similar size of state space is required for our method for exact calculation of the joint allele frequency spectrum in a structured coalescent with 2 islands.
<caption>**Table 1**</caption>

The six non-constant site patterns in an alignment of three labeled sequences ($S_1, S_2, S_3$). The alignment can be summarized in terms of the total number of singleton sites in the three sequences ($n_1, n_2, n_3$) and the number of shared mutations in the three sequences ($n_{12}, n_{13}, n_{23}$). We want to calculate the probability of a sample ($n_1, n_2, n_3, n_{12}, n_{13}, n_{23}$).

<table><thead><tr><th rowspan="2">Sequence</th><th colspan="6">Site patterns</th></tr><tr><th colspan="3">Singletons</th><th colspan="3">Doubletons</th></tr></thead><tbody><tr><td>$S_1$</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>$S_2$</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td></tr><tr><td>$S_3$</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td></tr><tr><td># sites</td><td>$n_1$</td><td>$n_2$</td><td>$n_3$</td><td>$n_{12}$</td><td>$n_{13}$</td><td>$n_{23}$</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 3. State transition diagram for a labeled tree with three sequences.</figcaption></figure>

## 5. More examples of the sampling formula

In this section we provide two examples of the sampling formula: a labeled sample with sample size $n = 3$ (Section 5.1) and a labeled two-island model with sample size $n = 2$ (Section 5.2). The examples allows for improved understanding of our population genetic sampling formula for the infinite sites mutation model and corresponding algorithm.

### 5.1. Labeled sample with sample size $n = 3$

The columns in an alignment of $n = 3$ sequences has 6 possible non-constant patterns as shown in Table 1. What are the key ingredients in our matrix-analytical method for calculating the probability of observing the patterns a number of times ($n_1, n_2, n_3, n_{12}, n_{13}, n_{23}$) in a sequence alignment? Here, $n_i$ is the number of singletons in sequence $i$, $i = 1, 2, 3$, and $n_{ij}$ is the number of shared mutations (doubletons) between sequence $i$ and $j$, $1 \le i < j \le 3$.

The ancestral process has four states (see Fig. 3). State 1 is the initial state with three singleton branches and zero doubleton branches. We enumerate this state (1, 1, 1, 0, 0, 0) where the first three entries indicate the presence of the three singleton branches, and the last three entries indicate the lack of doubleton branches. The chain jumps into State 2 if the next event is a coalescent between branch 1 and branch 2. After the merger the external branch 3 is still present, and a doubleton branch is present that could give rise to doubleton mutations in sequence 1 and sequence 2. We enumerate the state (0, 0, 1, 1, 0, 0). State 3 is enumerated (0, 1, 0, 0, 1, 0) and corresponds to a coalescent of branch 1 and branch 3, and state 4 is enumerated (1, 0, 0, 0, 0, 1) and corresponds to a coalescent between branch 2 and branch 3.

The total singleton branch length $B_i$ that can give rise to singletons in sequence $i$, $i = 1, 2, 3$, and the total doubleton branch length sequence $B_{ij}$ that can give rise to doubleton mutations in sequence $i$ and $j$, $1 \le i < j \le 3$, is multivariate phase-type distributed

$$B = (B_1, B_2, B_3, B_{12}, B_{13}, B_{23}) \sim \text{MPH}(\alpha, S, R)$$

with parameters

$$\alpha = (1, 0, 0, 0), \quad S = \begin{pmatrix} -3 & 1 & 1 & 1 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix} \quad \text{and} \quad R = \begin{pmatrix} 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}. \qquad (23)$$

Note that the reward matrix has rows that correspond to the enumeration of the states. We find

$$r_{\text{total}} = Re = (3, 2, 2, 2),$$

which corresponds to the total number of branches in each state. We get

$$P_{\text{total}} = (I - \lambda^{-1} \Delta (r_{\text{total}})^{-1} S)^{-1}$$
$$ = \begin{pmatrix} \frac{3\lambda}{(3+3\lambda)} & \frac{3}{(3+3\lambda)} & \frac{1}{3}(1+2\lambda) & \frac{2\lambda}{(3+3\lambda)} & \frac{3}{(3+3\lambda)} & \frac{1}{3}(1+2\lambda) & \frac{2\lambda}{(3+3\lambda)} & \frac{1}{3}(1+2\lambda) \\ 0 & \frac{2\lambda}{(1+2\lambda)} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{2\lambda}{(1+2\lambda)} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{2\lambda}{(1+2\lambda)} & 0 \end{pmatrix}, $$

and observe the nice interpretation of the various terms: in state 1 the two competing events are a mutation (rate $3\lambda$) and a coalescent (rate 3), each of the three trees are equally likely, and in state 2-4 the two competing events are a mutation (rate $2\lambda$) and a coalescent (rate 1). The exit probability vector is given by

$$ P_{\text{total}} = (I - P_{\text{total}})e = \left( \frac{3}{(3+3\lambda)} \frac{1}{(1+2\lambda)}, \frac{1}{(1+2\lambda)}, \frac{1}{(1+2\lambda)}, \frac{1}{(1+2\lambda)} \right), $$

and corresponds to the probability of being absorbed in the most recent common ancestor from the four states before any mutations occurs. The thinning vectors for the six mutation types are given by

$$ \begin{align*} r_1/r_{\text{total}} &= (1/3, 0, 0, 1/2), & r_2/r_{\text{total}} &= (1/3, 0, 1/2, 0), & r_3/r_{\text{total}} &= (1/3, 1/2, 0, 0), \\ r_{12}/r_{\text{total}} &= (0, 1/2, 0, 0), & r_{13}/r_{\text{total}} &= (0, 0, 1/2, 0), & r_{23}/r_{\text{total}} &= (0, 0, 0, 1/2), \end{align*} $$

and we get e.g.

$$ P_1 = P_{\text{total}} \Delta (r_1/r_{\text{total}}) = \begin{pmatrix} \frac{\lambda}{(3+3\lambda)} & 0 & 0 & \frac{1}{(3+3\lambda)} & \frac{\lambda}{(1+2\lambda)} \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{\lambda}{(1+2\lambda)} & 0 \end{pmatrix}. $$

The non-zero entries correspond to a singleton mutation in the external branch for sequence 1. The upper left entry is a mutation in state 1 given the starting state is state 1, the upper right entry is a mutation in state 4 given the starting state is state 1, and the lower right entry is a mutation in state 4 given the starting state is state 4. Another example of a thinned total sub-transition probability matrix is

$$ P_{12} = P_{\text{total}} \Delta (r_{12}/r_{\text{total}}) = \begin{pmatrix} 0 & \frac{1}{(3+3\lambda)} \frac{\lambda}{(1+2\lambda)} & 0 & 0 \\ 0 & \frac{\lambda}{(1+2\lambda)} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, $$

which corresponds to doubleton mutations in state 2 (i.e., the tree where branch 1 and branch 2 have coalesced).

It is interesting to observe that $P_{12}P_1 = 0$ because singleton mutations in sequence 1 are impossible after the merger of branch 1 and branch 2, but $P_1P_{12} \neq 0$ because the opposite is a possible event. We also have $P_{12}P_{13} = P_{12}P_{23} = P_{13}P_{23} = 0$ because the three types of doubletons are only possible in one of the trees.

## 5.2. Labeled two-island model with sample size $n = 2$

As a last example in this section, we consider the symmetric 2-island model with a sample size of $n = 2$. Let state 1 be the state where the alleles are on the same island, and state 2 be the state where the alleles are on different islands. Furthermore, let the migration rate from one island to the other be $m$. The joint branch length distribution $B = (B_1, B_2)$ of the left and right branch length, respectively, is then

$$ B \sim \text{MPH}(\alpha, S, R) $$

with parameters

$$ S = \begin{pmatrix} -2m-1 & 2m \\ 2m & -2m \end{pmatrix} \quad \text{and} \quad R = \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}, $$

where the reward matrix indicates that the branches accumulate mutations on the islands. We find

$$ r_{\text{total}} = Re = (2, 2), $$

which corresponds to the number of branches in each state. Using Section 4.1 we get

$$ P_{\text{total}} = \begin{pmatrix} 1 + (2m + 1)\frac{1}{2\lambda} & -\frac{m}{\lambda} \\ -\frac{m}{\lambda} & 1 + \frac{m}{\lambda} \end{pmatrix}^{-1} = \frac{1}{d(m, \lambda)} \begin{pmatrix} (2\lambda)(2m + 2\lambda) & (2m)(2\lambda) \\ (2m)(2\lambda) & (2\lambda)(2m + 1 + 2\lambda) \end{pmatrix}, $$

where $d(m, \lambda) = (2\lambda)^2 + (2\lambda) + 2(2\lambda)(2m) + (2m)$, and

$$ P_{\text{total}} = \frac{1}{d(m, \lambda)} (2m + 2\lambda, 2m). $$

We get the thinning vectors $r_1/r_{\text{total}} = r_2/r_{\text{total}} = (1, 1)/2$, and $P_1 = P_2 = \frac{1}{2}P_{\text{total}}$. Putting all this together gives the final sampling probability

$$ P(m_1, m_2) = \alpha P_{\text{total}}^{m_1+m_2} P_{\text{total}} \begin{pmatrix} m_1+m_2 \\ m_1 \end{pmatrix} \left(\frac{1}{2}\right)^{m_1+m_2}, \quad \text{where } (m_1, m_2) \in \mathbb{N}_{\ge 0}^2. $$
This formula reflects that given the number of mutations $m_1 + m_2$, the distribution of the mutations in the two branches is a symmetric binomial distribution.

In all previous examples, the initial state $\alpha$ has always been the first state. In the symmetric 2-island model we distinguish between the case where the two alleles are sampled from the same island ($\alpha_s = (1,0)$) or from different islands ($\alpha_d = (0,1)$). In the special case of zero mutations $m_1 = m_2 = 0$ and sampling from the same island we get

$$P_s = \alpha_s P_{\text{total}} = \frac{2m + 2\lambda}{d(m, \lambda)} = \frac{M + \theta}{\theta^2 + \theta + 2\theta M + M},$$

and from different islands we get the probability of identical alleles

$$P_d = \alpha_d P_{\text{total}} = \frac{2m}{d(m, \lambda)} = \frac{M}{\theta^2 + \theta + 2\theta M + M}.$$

Here, in both equations, we have used $\theta = 2\lambda$ and $M = 2m$ in the last equality because with this parametrization, the two equations are identical to equation (20) and equation (21) in Hudson (1990).

# 6. Application: Investigation of the site frequency spectrum

## 6.1. Sampling probabilities for the Kingman coalescent

If $m_{\text{total}}$ is the total number of mutations and $d$ is the number of mutation types ($d = n - 1$ for the SFS) we have that the number of possible configurations ($m_1, \dots, m_d$) equals all non-negative integer solutions to

$$m_1 + \dots + m_d = m_{\text{total}},$$

and this number equals the binomial coefficient $\binom{m_{\text{total}}+d-1}{m_{\text{total}}}$. For example, for configurations of the SFS with sample size $n = 6$ and $m_{\text{total}} = 10$, we have $d = n - 1 = 5$, and the number of possible configurations is equal to $\binom{10+4}{10} = \binom{14}{10} = 1001$. In our application we choose $n = 6$ and $m_{\text{total}} = 10$ because 1001 is a reasonable number of points to plot and configurations to consider. With sample size $n = 10$ and $m_{\text{total}} = 10$ we get $\binom{10+8}{10} = \binom{18}{10} = 43.758$ configurations, with $n = 6$ and $m_{\text{total}} = 20$ we get $\binom{20+4}{20} = \binom{24}{20} = 10.626$ configurations, and with $n = 10$ and $m_{\text{total}} = 20$ we get $\binom{20+8}{20} = \binom{28}{20} = 3.108.105$ configurations. The number of configurations thus increases very fast with sample size $n$ and total number of mutations $m_{\text{total}}$.

In Fig. 4 we summarize the sampling probabilities of the 1001 possible site frequency spectra with sample size $n = 6$ and total number of mutations $m_{\text{total}} = 10$. We condition on a sample size $n = 6$ and a total number of mutations $m_{\text{total}} = 10$, which makes Watterson's estimator $\hat{\theta}_W$ of $\theta$ ($\hat{\theta}_W = m_{\text{total}}/(\sum_{i=1}^{n-1} 1/i)$, see e.g. Section 4.1.1 in Wakeley (2008)) a constant. We sorted the compositions according to their probability (high to low) and normalized such that the probabilities sum to 1. In Fig. 4A we have plotted the cumulative probability of the compositions. It is perhaps expected but nevertheless striking that the spectra have huge variability in likelihood. The 10 most likely spectra account for more than 10% of the probability mass, the 60 most likely spectra account for more than 50% of the probability mass, and the 600 least likely spectra account for less than 1% of the probability mass.

The ten most likely spectra are summarized in Fig. 4B. As expected, they all contain many singletons and doubletons. Somewhat surprisingly, the four most likely site frequency spectra only contain singletons and pentagons. This means that the four most likely site frequency spectra have the smallest possible value of average pairwise distance between the $n = 6$ samples, namely $10 \cdot (1/3) = 3.33$ (each of the 10 mutations contribute $1/3$ to the average distance because the sample with the mutation is different from the other 5 samples and we have a total of $\binom{6}{2} = 15$ comparisons, and $5/15 = 1/3$).

In Fig. 4C we plot the joint average pairwise distance and sampling probability for each of the 1001 sample configuration. We have also indicated the number of singletons in the sample. The least likely spectra tend to have very few singletons. In Fig. 4D we have summarized the ten least likely spectra. They all contain zero singletons and are all extremely unlikely with probabilities less than $1.2 \cdot 10^{-10}$. We emphasize that evaluating such a small probability from simulations would require a huge number of replicates or a very clever sampling procedure. Our matrix-analytical sampling formula, which is highly precise and avoids simulation, is thus very useful in this context.

A classical statistic for investigating deviation from the standard $n$-coalescent is Tajima's $D$ (Tajima, 1989). Tajima's $D$ is defined by (e.g. equation (4.35) in Subsection 4.3.1 in Wakeley (2008))

$$D = \frac{\pi - m_{\text{total}}/\{\sum_{i=1}^{n-1} 1/i\}}{\text{sd}(n, m_{\text{total}})},$$

where $\pi$ is the average pairwise distance between the samples, $m_{\text{total}}$ is the total number of segregating sites, and $n$ is the sample size. The function $\text{sd}(n, m_{\text{total}})$ is the estimated standard deviation, and only depends on the sample size $n$ and the total number of segregating sites $m_{\text{total}}$. We observe that for a fixed sample size $n$ and number of segregating sites $m_{\text{total}}$ we have that Tajima's $D$ and $\pi$ are linearly dependent. In our analysis we have fixed $n = 6$ and $m_{\text{total}} = 10$, and therefore Tajima's $D$ and average pairwise distance differentiates between the samples in an identical manner. In Fig. 4C we compare the average pairwise distance and the sampling probability. It is clear from the plot that the average pairwise distance is a poor predictor of the probability of a composition. In particular, Fig. 4C suggests that reducing a sample to its Tajima's $D$ value may imply some strong loss of information when trying to reject Kingman's coalescent. Indeed, samples falling in the left top corner or right top corner of Fig. 4C have quite extreme Tajima's $D$ values and would thus lead it to reject this model, although they are among the most likely ones according to our sampling formula. A rejection test of the Kingman's coalescent based on our sampling formula could thus potentially outperform a test based on Tajima's $D$.

This suggestion would need to be explored more carefully, and the results would likely depend on the alternative model. In Fig. 5 we show the probability mass function of average pairwise distance. We observe that the distribution function is skewed to the right in accordance with Fig. 4B, C, and D. This result indicates that Tajima's $D$ is more powerful as a test for neutrality when the alternative model tends to produce large pairwise distances, but this would need further exploration because the detection power will depend on the demographic model, the parameters of the model, and on the sample size (Ferretti et al., 2010).
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 4. Conditional probabilities for all 1001 site frequency spectra with sample size *n* = 6 and total number of mutations *m*<sub>total</sub> = 10. A. Cumulative probability of all 1001 spectra.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>B. Ten most likely spectra. C. Joint values of average pairwise distance and probability of the composition. D. Ten least likely spectra.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 5. Probability distribution function for average pairwise distance with a fixed sample size *n* = 6 and number of mutations *m*<sub>total</sub> = 10.</figcaption></figure>

## 6.2. Sampling probabilities for the beta-coalescent

We also consider the Beta-coalescent. In particular we want to investigate how the sampling probabilities varies for different parametrizations of the Beta-coalescent. The rich class of multiple merger coalescent models was recently reviewed in Birkner and Blath (2022) (Chapter 8 in Baake and Wakolbinger (2022)). In this application we focus on the Beta-coalescent in the context of the site frequency spectrum.

In Table 2 we illustrate the rate matrix for the Beta-coalescent with sample size *n* = 6. Compared to the Kingman coalescent, the possibility of multiple mergers results in 11 new possible transitions between the transient states, and 7 new possible transitions directly to the MRCA.

We now describe the rate matrix for the block counting process of the Beta-coalescent with parameter *α* and sample size *n*. A state *a* = (*a*<sub>1</sub>, ..., *a*<sub>*n*−1</sub>) is encoded by all non-negative integer configurations where Σ<sub>*i*=1</sub><sup>*n*−1</sup> *ia*<sub>*i*</sub> = *n* with the interpretation that *a*<sub>*i*</sub> gives the number of lineages
<caption>**Table 2**</caption>

State space and illustration of the rate matrix for the block counting process with sample size $n = 6$. The binary mergers allowed in the Kingman's coalescent are marked with a '$k$', and additional mergers allowed in the Beta-coalescent are marked with a '$m$'. Empty entries are zero in the rate matrix.

<table><thead><tr><th rowspan="2">State</th><th rowspan="2">Number of<br>branches</th><th rowspan="2">State<br>number</th><th colspan="10">State number</th><th rowspan="2">MRCA</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td>(6,0,0,0,0)</td><td>6</td><td>1</td><td>.</td><td>k</td><td>m</td><td></td><td>m</td><td></td><td>m</td><td></td><td>m</td><td></td></tr><tr><td>(4,1,0,0,0)</td><td>5</td><td>2</td><td>.</td><td>k</td><td>k</td><td>k</td><td>m</td><td>m</td><td>m</td><td>m</td><td>m</td><td></td></tr><tr><td>(3,0,1,0,0)</td><td>4</td><td>3</td><td>.</td><td></td><td>k</td><td>k</td><td></td><td>m</td><td></td><td>m</td><td>m</td><td>m</td></tr><tr><td>(2,2,0,0,0)</td><td>4</td><td>4</td><td>.</td><td></td><td>k</td><td>k</td><td>k</td><td>m</td><td>m</td><td></td><td>m</td><td>m</td></tr><tr><td>(2,0,0,1,0)</td><td>3</td><td>5</td><td>.</td><td></td><td></td><td></td><td>k</td><td>k</td><td></td><td>m</td><td></td><td></td></tr><tr><td>(1,1,1,0,0)</td><td>3</td><td>6</td><td>.</td><td></td><td></td><td></td><td>k</td><td>k</td><td>k</td><td>m</td><td></td><td></td></tr><tr><td>(0,3,0,0,0)</td><td>3</td><td>7</td><td>.</td><td></td><td></td><td></td><td></td><td>k</td><td></td><td>m</td><td></td><td></td></tr><tr><td>(1,0,0,0,1)</td><td>2</td><td>8</td><td>.</td><td></td><td></td><td></td><td></td><td></td><td></td><td>k</td><td></td><td></td></tr><tr><td>(0,1,0,1,0)</td><td>2</td><td>9</td><td>.</td><td></td><td></td><td></td><td></td><td></td><td></td><td>k</td><td></td><td></td></tr><tr><td>(0,0,2,0,0)</td><td>2</td><td>10</td><td>.</td><td></td><td></td><td></td><td></td><td></td><td></td><td>k</td><td></td><td></td></tr></tbody></table>

subtending *i* leaves. Let *S* be the rate matrix between the states. It is possible to jump from a state *a* to a state *b* if exactly one new branch is created. Let *c* = *a* − *b*. Then a new branch is created if

$$c_{-} = - \sum_{i=1}^{n-1} c_i \mathbf{1}(c_i < 0) = 1, \qquad (24)$$

where $\mathbf{1}(\cdot)$ is the indicator function. For example, in Table 2 we can jump from state 3 (3,0,1,0,0) to state 8 (1,0,0,0,1) because $c = (2,0,1,0,-1)$.

The transition corresponds to a merger of two of the three singleton branches and the tripleton branch.

In the Beta-coalescent, the rate of merging $\ell$ out of $k$ branches is determined by a parameter $\alpha$ with support $1 \le \alpha < 2$, and is given by

$$\lambda_{k\ell} = \frac{\beta(\ell - \alpha, k - \ell + \alpha)}{\beta(\alpha, 2 - \alpha)}, \quad \ell = 2, \dots, k,$$

where $\beta(\cdot, \cdot)$ is the beta-function. For $\alpha \to 2$ we have $\lambda_{k\ell} \to 1$ for $\ell = 2$ and $\lambda_{k\ell} \to 0$ for $\ell = 3, \dots, k$, and we recover the Kingman coalescent. For $\alpha = 1$ we obtain the Bolthausen-Sznitman coalescent (Bolthausen and Sznitman, 1998).

We have that

$$c_{+} = \sum_{i=1}^{n-1} c_i \mathbf{1}(c_i > 0) \qquad (25)$$

of the branches merge. Furthermore, the number of possible ways to choose the branches in the multiple merger event is $\prod_{i=1}^{n-1} \binom{a_i}{c_i}$, where we define the binomial coefficient $\binom{a_i}{c_i}$ to be 1 if $j \le 0$ or $i = 0$. The final rate of coalescent from a state *a* to a state *b* is therefore

$$S_{a,b} = \mathbf{1}(c_{-} = 1) \lambda_{a_{+},c_{+}} \prod_{i=1}^{n-1} \binom{a_i}{c_i},$$

where $c = a - b$, $c_{-}$ is defined in (24), $c_{+}$ is defined in (25), and $a_{+} = \sum_{i=1}^{n-1} a_i$ is the number of branches in state *a*. For the Kingman coalescent $\lambda_{a_{+},c_{+}} = 1$ if $c_{+} = 2$ and $a_{+} \ge 2$, and zero otherwise.

In Fig. 6A we show the cumulative probability for the site frequency spectra for the Beta-coalescent with $\alpha \in \{1, 1.5, 2\}$, sample size $n = 4$ and total number of mutations $m_{total} = 10$. For small $\alpha$ we observe a more rapid increase of the cumulative probability of the sorted spectra than for large $\alpha$. This indicates that the Beta-coalescent with small $\alpha$ has a few compositions with large sampling probabilities compared to the Beta-coalescent with large $\alpha$, and this is indeed what we see in Fig. 6B where we show the sampling probabilities for compositions with a probability of at least 0.04.

# 7. Discussion

We have derived a new sampling formula for the probability of sampling population genetic data under a time-homogeneous coalescent process and the infinite sites mutation model. The sampling formula has an attractive interpretation in terms of a thinning of the total number of mutations. Furthermore, the sampling formula is easy to implement. In this final section we discuss future potential computational improvements and approximations, relaxation of the homogeneity assumption, and the possibility of using the sampling formula for improved simulation-based inference in population genetics.

## 7.1. Computational improvements and approximations

Our new sampling formula is numerically very stable and robust because the matrix products in Section 4.1 are products of (sub)transition probability matrices, and therefore all entries are positive and rows sum to at most one. Recall from Section 4 that the number of summands for a sample $\mathbf{m} = (m_1, \dots, m_d)$ in the sampling formula is given by the multinomial coefficient $c_\sigma = \binom{m_{total}}{m_1 \dots m_d}$, where $m_{total} = m_1 + \dots + m_d$ is the total number of mutations. This number quickly becomes very large, and this is the reason why the sampling formula is mainly useful in applications with a limited number of mutations. Recall from Section 4.1 that the number of matrix multiplications in our sampling formula is $c_\sigma(r_\sigma - 1)$, where $r_\sigma = m_{total}$ is the total number of mutations and $c_\sigma$ is the number of summands for the sample configuration. The number of summands can vary substantially between different configurations with the same sample size and total number of mutations.
<figure><img src="image_12.png" /><figcaption>Fig. 6. Conditional probabilities for the 66 site frequency spectra with sample size n = 4 and total number of mutations m<sub>total</sub> = 10. A. Cumulative probability of all 66 spectra for three different values of the free parameter α. B. Probabilities for the most likely spectra as a function of α. We show the probabilities for the 13 configurations with a probability larger than 0.04 for some value of α.</figcaption></figure>

Consider a configuration $(m_1, \dots, m_d)$ with $m_{total} = \sum_{j=1}^{d} m_j$ the total number of mutations. Then the number of summands (permutations) is given by the multinomial coefficient

$$c_{\sigma} = \binom{m_{\text{total}}}{m_1 \cdots m_d} = \frac{m_{\text{total}}!}{m_1! \cdots m_d!},$$

where $\sigma = \sigma(m_1, \dots, m_d)$. The multinomial coefficient is one if $m_i = m_{total}$ and $m_j = 0$ for $j \neq i$. If $m_{total} = kd$ then the largest value of the multinomial coefficient is obtained when $m_i = k$, $k = 1, \dots, d$. In this case we get

$$c_{\sigma(k, \dots, k)} = \frac{(kd)!}{(k!)^d} \qquad (26)$$

For the SFS with $n=6$ and $m_{total} = 10$ we get $d = n-1=5$ and $k=2$, and the multinomial coefficient is as large as $c_{\sigma(2,2,2,2,2)} = 113.400$. It thus takes more than 100.000 times longer to calculate the configuration (2, 2, 2, 2, 2) compared to the configuration (10, 0, 0, 0, 0)!

Using Stirling's formula, the largest multinomial coefficient (26) can be approximated by

$$c_{\sigma(k, \dots, k)} = \frac{(kd)!}{(k!)^d} \approx \frac{\sqrt{2\pi}\sqrt{k}\sqrt{d}d^{kd}}{\sqrt{2\pi}k^{kd/2}} = (2\pi k)^{(1-d)/2}d^{kd+1/2},$$

which reveals an exponential increase in k (i.e., the total number of mutations), and a faster than exponential increase in the mutation types d (i.e., in the sample size $n = d+1$). Finally, besides the number of matrix multiplications discussed above, the numerical complexity of our sampling formula is determined by the cost of each matrix multiplication, which is $p^3$, where p is the size of the state space. The state space size p is determined by the sample size n and possibly by certain characteristics of model (e.g. the number of demes in a structured model), and may also increase very fast for large samples, as discussed in Section 6.1.

There are, however, possibilities for speeding up and taking advantage of the structure of the sampling formula. Consider the problem of calculating the sampling probability of the most likely site frequency spectrum m = (7, 0, 0, 0, 3) for a sample size n = 6 and a total number of mutations m<sub>total</sub> = 10 (recall Fig. 4). Assume Step 1 and Step 2 in Section 4.2 have been carried out and we now want to calculate the sum of the matrices

$$\sum_{j=1}^{c_{\sigma}} \prod_{i=1}^{r_{\sigma}} \mathbf{P}_{\sigma_{ij}} = \mathbf{P}_1^7 \mathbf{P}_5^3 + \mathbf{P}_1^6 \mathbf{P}_5 \mathbf{P}_1 \mathbf{P}_5^2 + \mathbf{P}_1^6 \mathbf{P}_5^2 \mathbf{P}_1 \mathbf{P}_5 + \dots + \mathbf{P}_5^3 \mathbf{P}_1^7 \\ = (\mathbf{P}_1^2)^3(\mathbf{P}_1 \mathbf{P}_5)(\mathbf{P}_5^2) + (\mathbf{P}_1^2)^3(\mathbf{P}_5 \mathbf{P}_1)(\mathbf{P}_5^2) + (\mathbf{P}_1^2)^3(\mathbf{P}_5^2)(\mathbf{P}_1 \mathbf{P}_5) + \dots + (\mathbf{P}_5^2)^3(\mathbf{P}_5 \mathbf{P}_1)(\mathbf{P}_1^2)^3.$$

We observe (as indicated) that the sum is a function of the four matrix products $P_1^2, P_1P_5, P_5P_1$ and $P_5^2$. By pre-computing these matrix products we would only need half as many matrix multiplications in the sum of the matrices as before. This simple idea is actually even more powerful. For example recall from Section 3.4 that doubleton mutations are impossible after triplet mutations when the sample size is n = 4. In this case $P_3P_2 = 0$ and all the terms in the sum with this product can be removed.

Further strategies to speed up computations are to precompute products for ordered triplets, or even larger k-tuples of matrices $P_i$. This reduces the number of required matrix multiplications and k-tuples that have matrix products of zero can be removed. We also note that the sampling formula can easily be parallelized, as each of the summands can be computed independently of each other.

Finally, the order of the mutations also matters. Some orderings of mutations are unlikely, and an approximative calculation can be carried out where the terms corresponding to unlikely patterns of mutations are removed.

## 7.2. Relaxing the homogeneity assumption

It is natural to ask if the framework can be extended to time-inhomogeneous coalescent models. For example we might ask if the framework can be extended to a single population with a varying population size. Recall that our mathematical framework takes advantage of the analytical expression of the probability generating function (PGF) in Eq. (3). Unfortunately, to the best of our knowledge, an analytical expression for the
PGF is not available for time-inhomogeneous models. In fact, it is hard to imagine that it would be available because the reward transformation (scaling of exponential distributions) heavily depends on time-homogeneity.

A notable exception is demographic models with piecewise constant past population size variations where changes in coalescence rate exactly matches to transitions between ancestral coalescent states. Such model has proven to be useful in practice to infer past population sizes variations from the genome wide average SFS as in Liu and Fu (2015). However, even if this model could easily be implemented in our sampling formula, using the exact SFS at multiple independent non-recombining loci instead of the genome-wide expected SFS will lead to a set of locus-specific histories of past population size changes, and it is unclear how to reconcile these histories to obtain the multi-locus demographic history.

### 7.3. Improved Monte Carlo-based inference

Likelihood-based inference using coalescent simulations to estimate the likelihood of a genetic sample under various demographic models was initiated in a series of papers by Peter Beerli, Joseph Felsenstein and Mary Kuhner (see Chapter 27 in Felsenstein, 2004) and Robert Griffiths and Simon Tavare (see Tavaré, 2004). The latter approach uses importance sampling where the proposal distributions are derived from basic principles about the *n*-coalescent (Stephens and Donnelly, 2002; De Iorio and Griffiths, 2004a,b; De Iorio et al., 2005; Hobolth et al., 2008). The efficiency of the proposal determines the practical applicability of these algorithms, and they are often strongly limited by long computation times (Tewari and Spouge, 2015 and Rousset et al., 2018). One strategy to reduce computational time is to end the coalescent simulations before the MRCA when a state with two lineages is reached, and then use analytical probability of identity of pairs of genes for the last coalescent step (Rousset et al., 2018). Under the infinite sites mutation model, one could gain more computational time by stopping the coalescent simulation even earlier, i.e. with more than two lineages (perhaps even up to 6 ancestral lineages), and use our new sampling formulae to compute the probability of the remaining ancestral sample.

It could also be interesting to investigate if our new sampling formula for the infinite sites mutation model can be used to formulate improved sampling schemes for coalescent models. Such improved proposal distributions could also be applied for importance sampling beyond the classical *n*-coalescent framework and to the more general class of Lambda-coalescents (Birkner et al., 2011) or to the coalescent with recombination (Jenkins and Griffiths, 2011).

## CRediT authorship contribution statement

**Asger Hobolth:** Writing – original draft, Visualization, Software, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization. **Simon Boitard:** Writing – review & editing, Visualization, Funding acquisition, Conceptualization. **Andreas Futschik:** Writing – review & editing, Visualization, Funding acquisition. **Raphael Leblois:** Writing – review & editing, Visualization, Funding acquisition, Conceptualization.

## Funding information

We thank our institutions for supporting several research visits between the four of us. Two research visits of AH in Montpellier benefited from a travel grant funded by the LabEx CeMEB (ANR-10-LABX-04-01), Aarhus University Research Foundation, Denmark funded a visit of AF to visit AH in Aarhus, and JKU supported a visit of AH to AF in Linz.

AH is supported by the Novo Nordisk Foundation, Denmark grant number 22OC0079957. RL and SB were supported by the Occitanie Regional Council's program "Key challenge BiodivOc" managed by the University of Montpellier (DevOCGen project), and by recurrent funding from INRAE.

## Acknowledgments

We are grateful to Jere Koskela for fruitful discussions and Samuel Johnston for pointing out the relation of our sampling formula to the Faà de Bruno formula. We thank two anonymous reviewers and Ellen Baake for insightful and constructive comments that greatly helped improve the presentation.

## Data availability

No data was used for the research described in the article.

## References

* Baake, E., Wakolbinger, A., 2022. Probabilistic Structures in Evolution. EMS Press.

* Bhaskar, A., Wang, Y., Song, Y., 2015. Efficient inference of population size histories and locus-specific mutation rates from large-sample genomic variation data. Genome Res. 25 (2), 268–279.

* Birkner, M., Blath, J., 2022. Genealogies and inference for populations with highly skewed offspring distributions. pp. 151–158, Vol. Chapter 8.

* Birkner, M., Blath, J., Steinruecken, M., 2011. Importance sampling for Lambda-coalescents in the infinitely many sites model. Theor. Popul. Biol. 79, 155–173.

* Bisschop, G., 2022. Graph-based algorithms for Laplace transformed coalescence time distributions. PLoS Comput. Biol. 9 (18), e1010532.

* Bladt, M., Nielsen, B., 2017. Matrix-exponential distributions in applied probability. Matrix-exponential Distributions in Applied Probability, Springer.

* Bolthausen, E., Sznitman, A., 1998. On Ruelle's probability cascades and an abstract cavity method. Comm. Math. Phys. 197, 247–276.

* Chen, H., 2012. The joint allele frequency spectrum of multiple populations: a coalescent theory approach. Theor. Popul. Biol. 81 (2), 179–195.

* De Iorio, M., Griffiths, R.C., 2004a. Importance sampling on coalescent histories. I. Adv. in Appl. Probab. 36 (2), 417–433.

* De Iorio, M., Griffiths, R.C., 2004b. Importance sampling on coalescent histories. II: Subdivided population models. Adv. in Appl. Probab. 36 (2), 434–454.

* De Iorio, M., Griffiths, R.C., Leblois, R., Rousset, F., 2005. Stepwise mutation likelihood computation by sequential importance sampling in subdivided population models. Theor. Popul. Biol. 68 (1), 41–53.

* Ewens, W., 1972. The sampling theory of selectively neutral alleles. Theor. Popul. Biol. 3, 87–112.

* Excoffier, L., Dupanloup, I., Huerta-Sánchez, E., Sousa, V.C., Foll, M., 2013. Robust demographic inference from genomic and SNP data. PLoS Genet. 9 (10), e1003905.

* Felsenstein, J., 2004. Inferring Phylogenies. Sinauer Associates, Inc.

* Ferretti, L., Perez-Enciso, M., Ramos-Onsins, S., 2010. Optimal neutrality tests based on the frequency spectrum. Genet. 186 (1), 353–365.
Ganapathy, G., Uyenoyama, M.K., 2009. Site frequency spectra from genomic SNP surveys. Theor. Popul. Biol. 75 (4), 346–354.

Griffiths, R.C., 1982. The number of alleles and segregating sites in a sample from the infinite-alleles model. Adv. in Appl. Probab. 14 (2), 225–239.

Griffiths, R.C., Tavare, S., 1998. The age of a mutation in a general coalescent tree. Stoch. Models 14, 273–295.

Griffiths, R.C., Tavare, S., 2018. Ancestral inference from haplotypes and mutations. Theor. Popul. Biol. 122, 12–21.

Grimmett, G., Stirzaker, D., 1992. Probability and Random Processes. Oxford University Press.

Gutenkunst, R.N., Hernandez, R.D., Williamson, S.H., Bustamante, C.D., 2009. Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genet. 5 (10), e1000695.

Hankin, R.K.S., West, L.J., 2007. Set partitions in R. J. Stat. Softw. Code Snippets 23 (2), 1–12.

Hardy, M., 2006. Combinatorics of partial derivatives. Electron. J. Combin. 13 (1), R1.

Hobolth, A., Bladt, M., Andersen, L.N., 2021. Multivariate phase-type theory for the site frequency spectrum. J. Math. Biol. 83 (6–7), 63.

Hobolth, A., Rivas-Gonzàlez, I., Bladt, M., Futschik, A., 2024. Phase-type distributions in mathematical population genetics: An emerging framework. Theor. Popul. Biol. 157, 14–32.

Hobolth, A., Siri-Jegousse, A., Bladt, M., 2019. Phase-type distributions in population genetics. Theor. Popul. Biol. 127, 16–32.

Hobolth, A., Uyenoyama, M.K., Wiuf, C., 2008. Importance sampling for the infinite sites model. Stat. Appl. Genet. Mol. Biol. 7 (1).

Hudson, R.R., 1990. Gene genealogies and the coalescent process. Oxf. Surv. Evol. Biol. 7, 1–44.

Jenkins, P.A., Griffiths, R.C., 2011. Inference from samples of DNA sequences using a two-locus model. J. Comput. Biol. 18 (1), 109–127.

Kern, A., Hey, J., 2017. Exact calculation of the joint allele frequency spectrum for isolation with migration models. Genet. 207 (1), 241–253.

Kingman, J.F.C., 1982. The coalescent. Stochastic Process. Appl. 13 (3), 235–248.

Kulkarni, V.G., 1989. A new class of multivariate phase type distributions. Oper. Res. 37 (1), 151–158.

Kumagai, S., Uyenoyama, M.K., 2015. Genealogical histories in structured populations. Theor. Popul. Biol. 102, 3–15.

Liu, X., Fu, Y., 2015. Exploring population size changes using SNP frequency spectra. Nature Genet. 5, 555–559.

Lohse, K., Harrison, R., Barton, N., 2011. A general method for calculating likelihoods under the coalescent process. Genet. 189 (3), 977–987.

Mathematica, 2024. Wolfram research, inc.. URL https://www.wolfram.com/mathematica, Champaign, IL.

Neuts, M., 1981. Matrix-Geometric Solutions in Stochastic Models: an Algorithmic Approach. Johns Hopkins University Press.

Rivas-Gonzàlez, I., Andersen, L., Hobolth, A., 2023. PhaseTypeR: an R package for phase-type distributions in population genetics. J. Open Source Softw. 8, 50–54.

Rosenberg, N.A., 2020. Fifty years of Theoretical Population Biology. Theor. Popul. Biol. 133, 1–12, Fifty years of Theoretical Population Biology.

Rousset, F., Beeravolu, C.R., Leblois, R., 2018. Likelihood computation and inference of demographic and mutational parameters from population genetic data under coalescent approximations. J. Fr. Stat. Soc. 159 (3), 142–166.

Spence, J.P., Kamm, J.A., Song, Y.S., 2016. The site frequency spectrum for general coalescents. Genet. 202 (4), 1549–1561.

Stephens, M., Donnelly, P., 2002. Inference in molecular population genetics. J. R. Stat. Soc. Ser. B Stat. Methodol. 62 (4), 605–635.

Tajima, S., 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genet. 123 (3), 585–595.

Tavaré, S., 2004. Ancestral Inference in Population Genetics. Springer.

Tewari, S., Spouge, J., 2015. Coalescent: An open-science framework for importance sampling in coalescent theory. PeerJ 3, e1203.

Uyenoyama, M.K., Takebayashi, N., 2004. A simple method for computing exact probabilities of mutation numbers. Theor. Popul. Biol. 65 (3), 271–284.

Uyenoyama, M.K., Takebayashi, N., Kumagai, S., 2019. Inductive determination of allele frequency spectrum probabilities in structured populations. Theor. Popul. Biol. 129, 148–159.

Wakeley, J., 2008. Coalescent Theory: An Introduction. Roberts & Company Publishers.