Codex @codex  0

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 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.

For a neutral sample of size $n\ge 2$ under the infinite sites mutation model, the Watterson estimator is ${\hat{\theta }}_W=S/a_{n-1}$, where $a_{n-1}={\sum }_{i=1}^{n-1}1/i$. Since the coalescent has $ES=\theta a_{n-1}$, this estimator is unbiased. It matches the observed segregating-site count to its expected value, but need not maximize the exact segregating-site likelihood under the neutral coalescent.