Past exam of the mathematics course of the University of Cambridge / 2006 / iii / Paper 45 / 4

Use standard neutral coalescent time units, so each pair of ancestral lineages merges at rate 1. With $j$ lineages, there are $\left(\genfrac{}{}{0ex}{}{j}{2}\right)$ possible pairs; consequently the successive waiting times are independent withDuring this epoch there are $j$ branches, each of length $T_j$. Hence the total branch length of a neutral coalescent is $L={\sum }_{j=2}^{n}j{T}_j$, and linearity of expectation givesThe root branch above the sample's common ancestor is excluded: its mutations would be shared by every chromosome and would not be segregating sites in the sample.

In the infinite sites mutation model, mutations occur as independent Poisson processes at rate $\theta /2$ per lineage per unit of the coalescent time used in part (a). Conditional on the genealogy, superposition over branches gives a Poisson distribution with mean $\theta L/2$. Since this depends only on $L$, conditioning on total length alone givesEach mutation occurs at a new site and below the common ancestor, so it contributes exactly one segregating site. The law of iterated expectation therefore gives, with $a_{n-1}={\sum }_{i=1}^{n-1}1/i$,Here $\theta$ is fixed when taking the expectation. If it is also random under a proper prior with finite mean, the unconditional expectation is $a_{n-1}E\theta$.

Assume $n\ge 2$, $k\ge 0$, a proper prior density $\pi$ supported on nonnegative mutation rates, and independence of the mutation-rate prior from the neutral genealogy. PutBy Bayes' theorem, the joint posterior isThe normalizing constant is the prior predictive probability $P(S=k)$ and must be positive. For $k=0$, interpret the factor $(\theta L/2{)}^0$ as 1 even at zero. If $n=1$, there are no segregating sites: $S=0$ surely, so that case supplies no information about $\theta$.

Use Bayesian rejection sampling for a segregating-site count. Independently propose $\theta \sim \pi$ and $T_j\sim \mathrm{Exp}({\lambda }_j)$ for $j=2,\dots ,n$, and calculate $L={\sum }_{j}j{T}_j$. Draw an independent uniform $U$ on $(0,1)$ and accept the whole proposed vector exactly whenRepeat after a rejection. The right side is a probability mass function value of a Poisson distribution and is therefore between zero and one. On acceptance, the probability density function of a proposed vector is its prior density multiplied by this acceptance probability, divided by the probability of acceptance. Part (c) shows that this is precisely $f(\theta ,T\mid S=k)$. Independent proposals and uniforms produce independent posterior draws. Equivalently, one could simulate a count from a Poisson distribution with mean $\theta L/2$ and retain the proposal when that count equals $k$; the uniform version avoids simulating that extra count.

The acceptance probability is the prior average of the Poisson likelihood:Thus the expected number of proposals per accepted draw is $1/Z_k$. A numerical rate depends on the specified prior, sample size and observed count; it is not determined by $k$ alone.

For a more explicit one-dimensional expression, condition on $\theta$. The probability generating function, integrating each independent exponential epoch, isLet $p_k(\theta )$ be its coefficient of $z^k$. Then $Z_k=\int \pi (\theta )p_k(\theta )d\theta$. This identifies the acceptance rate with a fully specified marginal likelihood average.

The rejection rule can use a tighter envelope than 1. For $k>0$, differentiate $\log (e^{-x}{x}^k/k!)=-x+k\log x-\log (k!)$. Its derivative $-1+k/x$ changes from positive to negative at $x=k$, soKeep the same joint prior proposal but accept with probability $w_k(\theta ,T)/M_k$. Multiplication of the acceptance function by this constant leaves the accepted probability density function unchanged. The improved rate isFor $k>0$, $M_k<1$, so this is a genuine improvement; for large $k$, Stirling's formula gives $M_k\sim (2\pi k{)}^{-1/2}$. This is the best global constant envelope for an unrestricted positive prior proposal, because $L$ has support $(0,\infty )$ and $\theta L/2$ can approach $k$. For $k=0$ there is no improvement from this constant scaling. Further gains can come from proposals concentrating near $\theta L/2=k$, or integrating out the times first, but such changes require the appropriate proposal-density ratio and envelope rather than an uncorrected alteration of the sampling law.

The accepted mutation parameters have marginal probability density functionEstimate this probability density function with a suitably smoothed histogram or a boundary-aware kernel density estimator from the accepted draws. On a parameter region where $\pi (\theta )>0$, divide the estimate by the known prior density and maximize:The unknown $Z_k$ is constant in $\theta$ and cancels. Maximizing the posterior density itself gives a posterior mode, not generally a maximum-likelihood estimate; the two coincide for a flat prior on a region containing the maximum. A prior excluding candidate maximizers cannot recover them through its posterior output. Boundary behavior also matters: for $k=0$, the likelihood is decreasing in $\theta$, so the unrestricted maximum is at $\theta =0$.

The exact likelihood provides a useful independent check. Put $W=L/2$. Its Laplace transform is ${\prod }_{i=1}^{n-1}i/(i+s)$, since $j{T}_j/2$ is exponential of rate $j-1$. The probability density functionhas that same transform: substitute $u=e^{-w}$ to obtain $(n-1){\int }_0^1{u}^s(1-u{)}^{n-2}du={\prod }_{i=1}^{n-1}i/(i+s)$. Thus, expanding its binomial factor and integrating each term,This is the segregating-site likelihood under the neutral coalescent. It can be maximized numerically to check the posterior-density estimate; the generating-function representation avoids cancellation in large alternating sums. For $n=2$, it reduces to ${\theta }^k/(1+\theta {)}^{k+1}$, whose maximum is at $\theta =k$ for $k>0$.

Alternatively, differentiating the integrated likelihood givesAn interior maximum therefore satisfies $\theta =k/E[W\mid S=k,\theta ]$. Conditional expectations of $W$ can be estimated by smoothing the accepted pairs $(\theta ,W)$ near each candidate parameter value, yielding another simulation-based score equation. The simple estimator $k/a_{n-1}$ is the unbiased Watterson estimator from part (b), and is not generally the likelihood maximizer.