Inferring poisson rate from interval determined by data

Question

I have a dataset of arrivals, which are from a Poisson process. For the purposes of this question, let's say they're arrivals of cars on a particular road. My goal is to infer the gamma posterior for the poisson rate parameter, $\lambda$.

However, the dataset was sampled peculiarly:

1. Events were counted until exactly $N$ had been measured
2. Only the absolute times ($t_i, ... t_N$) of events were recorded

I believe that it's incorrect to estimate $\lambda$ using $t_N - t_i$ as the length of the sample interval, because the start and end of the interval use knowledge of the times of datapoints. Is it actually incorrect? If so, why? (in terms of probability)

Would it instead be correct to estimate $\lambda$ from the inter-event-arrival times? that is, if $\tau = t_N - t_0$, then the posterior parameters (using the shape and rate formulation as described on wikipedia) are:

$\alpha' = \alpha + N - 1$

$\beta' = \beta + \tau$

Where $\alpha'$ and $\beta'$ are the updated parameters after observing evidence.

Answer 1

I think you are correct, with minor quibbles. The observed times must have been $(t_0, t_1, \dots, t_N).$ Your $N$ observed exponential interarrival times are $d_i = t_{i} - t_{i-1}$, for $i = 1, \dots, N$. Their sum is your $\tau$.

The likelihood function based on exponential interarrival times is $$\prod_{i=1}^N \lambda \exp(-\lambda d_i) = \lambda^N e^{-\lambda\tau}.$$ Thus your posterior is gamma with the parameters you give in your question.

As a reality check, I tried generating $N = 50$ observations from an exponential distribution with rate $\lambda = 1/5 = .2,$ and then cutting 2.5% from each tail of your posterior based on an almost noninformative prior GAMMA(shape = 1/2, rate = 1/100) and the data. With the simulated data from one run I got (0.171, 0.297) as a Bayesian probability interval. Several other runs gave intervals that covered the known $\lambda = .2.$ In case you are familiar with R, I'm pasting the code below.

lam = 1/5;  n = 50;  d = rexp(n, lam)
pr.sh = .5;  pr.rt = 1/100
qgamma(c(.025,.975), pr.sh + n, pr.rt + sum(d))
## 0.1705835 0.2969390

KW: Gamma-Gamma Conjugates, Bayesian Inference