= Solution
The accepted values $\theta^{(1)},\ldots,\theta^{(R)}$ have marginal <posterior density>
$$
h(\theta\mid k)=\int f(\theta,t\mid k)\,dt=\frac{\pi(\theta)\ell_k(\theta)}{m_k}.
$$
Estimate this one-dimensional density by a <histogram> or <kernel density estimation> from the simulated values. Where $\pi(\theta)>0$, <likelihood reconstruction from posterior simulation> gives
$$
\boxed{\widehat\theta_{\mathrm{ML}}\approx\underset{\theta}{\operatorname{argmax}}\ \frac{\widehat h_R(\theta\mid k)}{\pi(\theta)}.}
$$
Use 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:
$$
\frac{R^{-1}\#\{r:\theta^{(r)}\in I_b\}}{\int_{I_b}\pi(u)\,du}\ \longrightarrow\ \frac{\int_{I_b}\pi(u)\ell_k(u)\,du}{m_k\int_{I_b}\pi(u)\,du}.
$$
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. Put
$$
r_k(t)=\int_0^\infty\pi(u)a_k(u,t)\,du.
$$
Their <posterior density> is $q(t)r_k(t)/m_k$. Therefore, for any trial value $\vartheta$,
$$
\frac1R\sum_{r=1}^R\frac{a_k(\vartheta,T^{(r)})}{r_k(T^{(r)})}\ \longrightarrow\ \frac1{m_k}\int q(t)a_k(\vartheta,t)\,dt=\frac{\ell_k(\vartheta)}{m_k}.
$$
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.
Back to article page
