Solution (source code)

= Solution

Use the joint <prior> density $g(\theta,t)=\pi(\theta)q(t)$ as the proposal for <rejection sampling>. Since the Poisson <probability> $a_k(\theta,t)$ lies in $[0,1]$, the following algorithm is valid:

* Draw $\theta$ from the proper <prior distribution> $\pi$.
* Independently draw $T_j$ from the <exponential distributions> with rates $\binom j2$, for $j=2,\ldots,n$, and calculate $L=\sum jT_j$.
* Draw $U$ uniformly on $[0,1]$, independently of all previous draws. Retain $(\theta,T)$ if $U\le e^{-\theta L/2}(\theta L/2)^k/k!$; otherwise restart.

The joint density of a proposed pair together with its acceptance event is $g(\theta,t)a_k(\theta,t)$. Its integral is $m_k$, so conditioning on acceptance gives precisely $g a_k/m_k$, the desired <posterior density>. Repeating independent proposal trials until each acceptance produces independent <posterior> observations. This is <Bayesian rejection sampling for a segregating-site count>.

An equivalent implementation simulates $S^*\sim\operatorname{Poisson}(\theta L/2)$ and accepts exactly when $S^*=k$. Its acceptance <probability> is the same Poisson mass. \b[Propose from the joint <prior> and accept according to the exact observed-count <likelihood>.] No approximation to the <posterior distribution> is involved in the accepted draws.