Past exam of the mathematics course of the University of Cambridge / 2005 / iii / Paper 45 / 3

Use standard diploid coalescent time units of $2N_e$ generations, where $N_e$ is the fixed effective population size, and take $n\ge 2$. In the neutral Kingman's coalescent, each unordered pair of ancestral lineages merges at rate one. With $j$ lineages, the next coalescence time therefore has an exponential distribution with rateThe memoryless construction and the Markov property make these epoch durations independent. During epoch $j$ there are $j$ branches, each of duration $T_j$, so the total branch length of a neutral coalescent isBy linearity of expectation,The ancestral branch above the sample's most recent common ancestor is excluded: mutations there would be shared by all sampled chromosomes and would not create segregating sites.

Let $\mu$ be the mutation rate per chromosomal segment per generation, and use the population-scaled parameter $\theta =4N_e\mu$. In $2N_e$-generation coalescent units, mutations occur on each ancestral lineage at rate $\theta /2$. Independent mutations on all branches of a given phylogenetic tree superpose to give a Poisson distribution with mean $\theta L/2$.

Under the infinite sites mutation model, each mutation occurs at a distinct position, and every mutation below the sample's common ancestor separates a nonempty proper subset of sampled chromosomes from the remainder. Thus the mutation count is exactly the segregating site count, givingBy the law of total expectation and the previous branch-length formula,This also gives the Watterson estimator $S/a_n$ as an unbiased moment estimator of $\theta$. Here all sample segregating sites are counted, including singletons; filtering sites by a minimum sample frequency or by external SNP discovery would introduce genetic-study ascertainment and change this calculation.

Let $t=(t_2,\dots ,t_n)$, with all $t_j>0$, and put $L(t)={\sum }_{j=2}^{n}j{t}_j$. The independent exponential epoch times have joint densityTake a proper prior distribution density $\pi (\theta )$ on $\theta \ge 0$, independent of the neutral genealogy. Independence is appropriate because the specified neutral Kingman's coalescent does not depend on the mutation parameter. The conditional observed-count likelihood function isIts prior predictive probability, or Bayesian model evidence, isAssuming $m_k>0$, Bayes' theorem gives the joint posterior of mutation rate and coalescent times:It is zero outside the prior support and positive-time region. For $k=0$, use the usual convention $0^0=1$ in the Poisson mass. An improper prior cannot be sampled as the proposal in the following algorithm, even in cases where its formal posterior can be normalized.

Use the joint prior density $g(\theta ,t)=\pi (\theta )q(t)$ as the proposal for rejection sampling. Since the Poisson probability $a_k(\theta ,t)$ lies in $[0,1]$, the following algorithm is valid:

The joint density of a proposed pair together with its acceptance event is $g(\theta ,t)a_k(\theta ,t)$. Its integral is $m_k$, so conditioning on acceptance gives precisely $g{a}_k/m_k$, the desired posterior density. Repeating independent proposal trials until each acceptance produces independent posterior observations. This is Bayesian rejection sampling for a segregating-site count.

An equivalent implementation simulates $S^{\ast }\sim \mathrm{Poisson}(\theta L/2)$ and accepts exactly when $S^{\ast }=k$. Its acceptance probability is the same Poisson mass. Propose from the joint prior and accept according to the exact observed-count likelihood. No approximation to the posterior distribution is involved in the accepted draws.

Integrating the acceptance probability over the joint-prior proposal givesThe mean number of independent proposal trials per acceptance is $1/m_k$, by the geometric distribution. To make the acceptance rate explicit for an arbitrary specified prior, define the segregating-site likelihood under the neutral coalescent by ${\ell }_k(\theta )=E_T[a_k(\theta ,T)]$. ThenFor the last identity, $W=L/2$ is the sum of independent exponential variables with rates $1,\dots ,n-1$; combine their Laplace transforms with the conditional Poisson probability generating function. If an explicit mass is wanted, $W$ also has density $(n-1)e^{-w}(1-e^{-w}{)}^{n-2}$: it is the maximum of $n-1$ independent unit exponential variables, whose successive spacings have exactly these exponential rates. Expanding $(1-e^{-w}{)}^{n-2}$ and integrating the Poisson mass givesUse ${\theta }^0=1$ for $k=0$. The positive integral or generating-function expression may be numerically preferable to this alternating sum.

A sharper rejection envelope improves the algorithm without changing its proposal. For $k>0$, maximize the Poisson mass over $\lambda =\theta L/2$:It increases until $\lambda =k$ and decreases afterwards, so its maximum is $M_k=e^{-k}{k}^k/k!$. For $k=0$, the maximum is $M_0=1$. Replace the acceptance rule by $U\le a_k(\theta ,T)/M_k$. This is bounded by one, and its accepted density is unchanged because the factor $1/M_k$ cancels on normalization. ThusFor $k>0$ it is strictly more efficient, with improvement factor $1/M_k\sim \sqrt{2\pi k}$ by Stirling's formula. For $k=0$ this envelope improvement gives no increase. Evaluating the Poisson mass and using a uniform draw also avoids generating a Poisson count only to discard it. A proposal concentrated near the posterior can improve efficiency further, but then its density and a valid rejection bound must be included.

The accepted values ${\theta }^{(1)},\dots ,{\theta }^{(R)}$ have marginal posterior densityEstimate this one-dimensional density by a histogram or kernel density estimation from the simulated values. Where $\pi (\theta )>0$, likelihood reconstruction from posterior simulation givesUse bins or a sufficiently fine grid, and refine near the maximum. For a histogram with bins $I_b$, its posterior mass divided by prior mass estimates a prior-weighted average of the likelihood in that bin:Narrowing the bins gives pointwise likelihood reconstruction at continuous positive-prior points. A flat proper prior over the search interval makes the highest posterior-density bin approximately the highest-likelihood bin. With a nonconstant prior, maximizing the posterior alone gives a maximum a posteriori estimate, generally different from maximum likelihood estimation. Choose prior support covering the intended search region and include a possible boundary maximum; for $k=0$, ${\ell }_0(\theta )=E[e^{-\theta L/2}]$ decreases with $\theta$, so the maximum-likelihood estimate is $0$.

The sampled genealogies provide a second likelihood reconstruction when the one-dimensional prior integral is computable. PutTheir posterior density is $q(t)r_k(t)/m_k$. Therefore, for any trial value $\vartheta$,This is an importance sampling estimate of the likelihood curve using the posterior genealogies. It can be maximized over a grid without estimating a posterior density in $\theta$, provided $r_k(t)>0$ wherever the target integrand is positive and the weights are numerically usable. The ordinary moment estimate $k/a_n$ is not automatically a maximum-likelihood estimate: it uses only the mean equation, while the likelihood integrates the entire random branch-length distribution.