Solution (source code)

= Solution

Assume $n\ge2$, $k\ge0$, a proper <prior density> $\pi$ supported on nonnegative mutation rates, and <independence> of the mutation-rate prior from the neutral genealogy. Put
$$
q(t)=\prod_{j=2}^n\lambda_j e^{-\lambda_jt_j}\mathbf1_{\{t_j>0\}},\qquad
L(t)=\sum_{j=2}^n jt_j,\qquad
w_k(\theta,t)=e^{-\theta L(t)/2}\frac{(\theta L(t)/2)^k}{k!}.
$$
By <Bayes' theorem>, the joint posterior is
$$
\boxed{f(\theta,t\mid S=k)=\frac{\pi(\theta)q(t)w_k(\theta,t)}{Z_k},\quad
Z_k=\int_0^\infty\!\pi(u)\int_{(0,\infty)^{n-1}}q(t)w_k(u,t)\,dt\,du.}
$$
The 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$.