Reversible-jump MCMC and Poisson processes

Suppose we have a time interval $$t \in [0, T]$$ in which events occur as a Poisson process with some arbitrary time-dependent rate $$\lambda(t)$$. These events occur at times $$Y=(Y_1, Y_2, \dotso, Y_M)$$ for $$M$$ events, where $$Y$$ and $$M$$ are both random variables, and $$0 \leq Y_i \leq T$$. $$T$$ is fixed.

The objective is to estimate the number of events $$M$$ and their times $$Y$$ using MCMC.

We know that $$M$$ follows a Poisson distribution with rate $$g = \int\limits^T_0 \lambda(s) ds$$, and hence $$p(M|T) = \frac{e^{-g} g^M}{M!}$$.

I believe the joint probability density is:

$$p(Y,M|T) = p(Y|M,T) p(M|T) = \Big( \prod\limits_{i=1}^M \frac{\lambda(y_i)}{g} \Big) \Big( \frac{e^{-g} g^M}{M!} \Big) = \frac{e^{-g}}{M!} \prod\limits_{i=1}^M \lambda(y_i)$$

I have confirmed that the integral under this sums to 1, with respect to $$Y$$ and $$M$$.

However, I am unable to sample from this distribution using MCMC. Currently, I have two MCMC proposals, and these are selected with equal probabilities at each step in the chain.

1. Move event. Sample one of the indices $$i$$ in $$Y$$ and sample $$Y_i^\prime \sim \text{U}(0,T)$$. The Hastings ratio $$\frac{q(M,Y|M^\prime,Y^\prime)}{q(M^\prime,Y^\prime|M,Y)} = \frac{\frac{1}{MT}}{\frac{1}{MT}} = 1$$. If $$M=0$$ the proposal is rejected.

2. Create/destroy event.

a. With 0.5 probability, create an event. $$M^\prime \leftarrow M+1$$. Sample an index $$i \sim \text{U}(1, \dotso, M)$$ and shift all of the downstream indices up by one: $$Y_{i+1}^\prime \leftarrow Y_i, Y_{i+2}^\prime \leftarrow Y_{i+1}, \dotso$$. The new value is sampled uniformly $$Y_i^\prime \sim \text{U}(0, T)$$.

b. With 0.5 probability, destroy an event. $$M^\prime \leftarrow M-1$$. Sample an index $$i \sim \text{U}(1, \dotso , M)$$ and remove it, shifting all of the downstream indices down by one: $$Y_{i}^\prime \leftarrow Y_{i+1}, Y_{i+1}^\prime \leftarrow Y_{i+2}, \dotso$$. If $$M=0$$ the proposal is rejected.

The Hastings-ratio of a birth event is $$\frac{q(M,Y|M^\prime,Y^\prime)}{q(M^\prime,Y^\prime|M,Y)}=\frac{\frac{1}{2}\frac{1}{M+1}}{\frac{1}{2}\frac{1}{T}} = \frac{T}{M+1}$$, and $$\frac{M}{T}$$ for a death event.

But, as mentioned, this algorithm does not sample the correct distribution (i.e., the estimated number of events $$M$$ is generally larger than the expected value $$g$$). I believe the issue comes down to a Jacobian term used in the second operator, which involves a change in dimension. The acceptance probability during MCMC is thus:

$$\alpha = \min\Big(1, \frac{p(Y^\prime, M^\prime|T) q(M,Y|M^\prime,Y^\prime)}{p(Y, M|T) q(M^\prime,Y^\prime|M,Y)} |J| \Big)$$

for some Jacobian $$J$$. But what is $$J$$? Or perhaps the problem lies elsewhere. Any help appreciated.

• Welcome to CV, Jordan. I'm sorry about the unexplained downvote: it's clear to me you have put effort and research into your question (+1).
– whuber
Commented Jan 31 at 22:38
• Thanks! Happy to clarify if I missed any important details Commented Jan 31 at 22:52

First, a quick correction: step 2a should read "Sample an index $$i \sim U(1, 2, \dotso, M, M+1)$$", rather than $$M$$, because there are $$M+1$$ possible places a new element can be inserted in a list with $$M$$ elements.
$$\frac{q(M,Y|M',Y')}{q(M',Y'|M,Y)} = \frac{\frac{1}{2} \frac{1}{M+1}}{\frac{1}{2} \frac{1}{M+1} \frac{1}{T}} = T \text{ for a birth and } \frac{1}{T} \text{ for a death}.$$