Past exam of the mathematics course of the University of Cambridge / 2004 / iii / Paper 39 / 1

During the epoch with $j$ ancestral lineages, Kingman's coalescent has waiting time $T_j\sim \mathrm{Exp}((\genfrac{}{}{0ex}{}{j}{2}))$. In the stated time units, each lineage accumulates mutations at rate $Nu=\theta /2$, so the mutation count on a neutral coalescent tree decomposes asThe law of total expectation and the mean of an exponential distribution give $\mathbb{E}{Y}_j=(j\theta /2)/\left(\genfrac{}{}{0ex}{}{j}{2}\right)=\theta /(j-1)$. ThereforeThis calculation counts mutation events. Equating the count with observed segregating sites additionally uses the infinite sites mutation model.

Apply the law of total variance to each epoch. The Poisson distribution has equal conditional mean and variance, while $\mathrm{Var}(T_j)={\left(\genfrac{}{}{0ex}{}{j}{2}\right)}^{-2}$. HenceThe independent mutation counts in coalescent epochs permit variance additivity for independent random variables. ThusThe first term is conditional Poisson distribution variation; the second comes from the random lengths of the ancestral tree. These are the moments of the mutation count on a neutral coalescent tree.

The Watterson estimator isThe expectation above makes it an unbiased estimator. Its variance, given by the variance of the Watterson estimator, isIndeed, the harmonic sum satisfies $H_{n-1}\sim \log n$ and $H_{n-1}^{(2)}$ stays bounded. This proves mean-square convergence; the Chebyshev inequality then shows that it is a consistent estimator, since $P(|\hat{\theta }-\theta |>\epsilon )\le \mathrm{Var}(\hat{\theta })/{\epsilon }^2\to 0$.

Conditioning on the epoch length gives the probability generating functionThe Laplace transform of an exponential distribution then yieldsThis is the geometric mutation count in a coalescent epoch. With the zero-based geometric distribution convention,At $\theta =0$, this reduces to $Y_j=0$ almost surely.

For fixed $\theta >0$, use the Lyapunov condition for the independent centered epoch counts $Z_i=Y_{i+1}-\theta /i$. Their total variance is $s_n^2=\theta H_{n-1}+{\theta }^2{H}_{n-1}^{(2)}\sim \theta \log n$.

Write ${\mu }_i=\theta /i$. Differentiating the geometric distribution probability generating function gives $\mathbb{E}[Y_{i+1}^3]={\mu }_i+6{\mu }_i^2+6{\mu }_i^3$. Since $|Y_{i+1}-{\mu }_i{|}^3\le 4Y_{i+1}^3+4{\mu }_i^3$, the sum of the centered absolute third moments is $O(\log n)$. ConsequentlyThe Lindeberg-Feller central limit theorem applies through the Lyapunov condition, giving $(S-\theta H_{n-1})/s_n\stackrel{d}{\to }N(0,1)$. Now $\sqrt{\log n}{s}_n/H_{n-1}\to \sqrt{\theta }$, so the Slutsky theorem proves the asymptotic normality of the Watterson estimator:The limiting variance is $\theta$. When $\theta =0$, every count vanishes and the limiting law is the point mass at zero.

The asymptotic normality of the Watterson estimator gives approximate standard error $\sqrt{\theta /\log n}$. Thus precision at a fixed nonrecombining locus improves only logarithmically with sample size. To halve the asymptotic standard error, one needs $\log n^{\prime }\simeq 4\log n$, or roughly $n^{\prime }=n^4$, rather than four times as many sampled chromosomes.

For large samples and positive $\hat{\theta }$, a plug-in normal approximation gives $\hat{\theta }\pm 1.96\sqrt{\hat{\theta }/\log n}$ as an approximate 95% confidence interval. At low mutation counts it can perform poorly or give a negative endpoint; the exact count model is then preferable. Multiple independent, comparable loci can provide additional independent information, whereas tightly linked chromosomes still share the random ancestry at one locus. Genetic recombination, nonneutral evolution and changing population size require corresponding changes to the model.