2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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). $\endgroup$
    – whuber
    Commented Jan 31 at 22:38
  • $\begingroup$ Thanks! Happy to clarify if I missed any important details $\endgroup$
    – Jordan
    Commented Jan 31 at 22:52

1 Answer 1

0
$\begingroup$

I believe I have a solution now. There were two issues.

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.

Second, the creation proposal kernel was not taking the random indexing of this new event into account, but the reverse proposal (destroy) was. The corrected Hastings ratio of the birth-death proposal is:

$$ \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}. $$

Expectation of M and estimated number of events are now the same. Figure below is from 100 MCMC replicates across 100 expected values.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.