Bayesian parameter estimation of a Poisson process with change/no-change observations at irregular intervals

Question

Consider a Poisson process with unknown parameter $\lambda$.
We perform a sequence of $n$ observations at intervals $\overline{t}=t_1,\,t_2,\,\dots,\,t_n$. Each observation is a binary variable $x_i$ equal to zero if no changes occurred during interval $t_i$ or equal to one if one or more changes occurred:
$\forall i,$
$\Pr[X_i = 0 | \lambda,\,t_i] = \exp(-\lambda t_i)$
$\Pr[X_i = 1 | \lambda,\,t_i] = 1 - \exp(-\lambda t_i)$

The likelihood function for all observations is:
$P[\overline{x}| \lambda,\overline{t}] = (\prod_{x_i=0}\exp(-\lambda t_i)) \cdot (\prod_{x_i=1}(1-\exp(-\lambda t_i)))$

Assuming that the observation intervals are chosen independently from anything else, I would like to estimate the parameter $\lambda$ under some reasonable prior
$\lambda^\star = argmax_\lambda P[\lambda | \overline{x},\,\overline{t}]$
$P[\lambda | \overline{x},\,\overline{t}] \propto P[\overline{x}| \lambda,\overline{t}] \cdot P[\lambda]$

If possible, I would like to use a conjugate prior in order to perform recursive estimation, but I'm not sure whether one exists.

My questions:
1. Does a conjugate prior for that likelihood exist?
2. If a conjugate prior does not exist, what estimator can I use? I'm interested in an estimator that can be updated incrementally for each observation, without keeping track of all the observation history.

Answer 1

Notation

First I'll re-express the log-likelihood:

$$\log L(\lambda)=\sum_i \bigg[ x_i \log(1-e^{- \lambda t_i}) - \lambda (1-x_i)t_i\bigg].$$

Let

• $\tau_0 = \sum_i (1-x_i)t_i$,
• $\tau_1 = \sum_i x_i t_i$,
• $\tau_2 = \sum_i x_i t_i^2$,
• $n_1 = \sum_i x_i$.

Approximation

According to Wolfram Alpha,

$$ \log(1-e^{-x}) \approx \log(x) - \frac{x}{2} + \frac{x^2}{24}.$$

As $x \searrow 0$, $\log(1-e^{-x}) \rightarrow -\infty$, and as $x \rightarrow \infty$, $ \log(1-e^{-x}) \nearrow 0$. On the other hand, the approximation contains an upward-facing parabola, so it will diverge to $+\infty$ as $x \rightarrow \infty$. So when we find the maximum of the log-likelihood, only the smaller root of the parabola will be in the region where the approximation is accurate.

Estimation

Up to constant terms, substituting in the above approximation results in the following approximate log-likelihood:

$$\log \tilde{L}(\lambda) = n_1 \log(\lambda) - \lambda \bigg(\tau_0 + \frac{\tau_1}{2}\bigg) + \lambda^2\bigg(\frac{\tau_2}{24}\bigg).$$

The approximate MLE, $\hat{\lambda}$, satisfies

$$ \begin{array}{rcl}\frac{{\rm {d}}}{{\rm {d}\lambda}}\log\tilde{L}(\lambda)\bigg|_{\lambda=\hat{\lambda}} & = & 0\\n_{1}-\hat{\lambda}\bigg(\tau_{0}+\frac{\tau_{1}}{2}\bigg)+\hat{\lambda}^{2}\bigg(\frac{\tau_{2}}{12}\bigg) & = & 0\\ \Rightarrow \hat{\lambda} & = & \frac{6\tau_{0}}{\tau_{2}}+\frac{3\tau_{1}}{\tau_{2}}-\frac{6}{\tau_{2}}\sqrt{\bigg(\tau_{0}+\frac{\tau_{1}}{2}\bigg)^{2}-\frac{n_{1}\tau_{2}}{3}}.\end{array} $$

Prior distribution

The gamma distribution isn't exactly conjugate but will certainly be easy to combine with the above approximate log-likelihood.

Answer 2

There seems to be no conjugate prior for the Poisson regression model (see). However, you can use any proper prior on $\lambda$ that reflects your prior knowledge on it. For example: Gamma distribution, log-normal distribution, exponential distribution, among others. If your prior information is vague, then you can use a heavy-tailed prior, for instance a half-Cauchy distribution with scale parameter $\sigma=5$ (you may want to try with larger values).

In order to simulate from the posterior, you can employ the R package mcmc and the command metrop() where you simply have to input the log-posterior, up to a proportionality constant, and adjust the variance of the proposal to obtain a reasonable acceptance rate. If you want more details, please post your sample or a simulated sample.

Using the posterior sample you can calculate and compare several estimators such as Maximum a posteriori (MAP), the posterior mean, the posterior median, and etcetera. Although the use of these estimators may not allow for an explicit expression of the posterior, the simulation step is reasonably fast and therefore you can run a simulation everytime you add a new observation.

Answer 3

Not bayesian estimation

When we suppose that all yours intervals $t_i$ have equal size $\Delta$, we can use the following asymptotically unbiased estimator proposed in another stack.exchange question:

$$ \newcommand{\E}{\mathbb{E}} \hat\lambda = - \Delta^{-1} \log \frac{W}{N} $$

$$ \begin{aligned} \Delta & \text{ -- size of the time interval between observations} \\ W & \text{ -- number of intervals without changes} \\ N & \text{ -- total number of intervals} \end{aligned} $$ also we calculated its bias: $$ \E \left[ \hat\lambda - \lambda \right] \approx \frac{1}{2N\Delta} \left(e^{\lambda\Delta} - 1\right) $$