Estimating event duration by sampling with a Poisson process

Question

Take an event to be a set $E \subseteq {\mathbb R}^{\ge 0}$. We want to estimate the duration of the event (the measure of $E$) by sampling using a Poisson process $N$. That is, for some experiment duration $D$, at each time $t\leq D$ such that $\forall t'>t\colon N(t')>N(t)$, we check whether $t\in E$. See TagTime: Stochastic Time Tracking for Space Cadets for more context.

Given $n$, the number of sampled times $t\in E$, and $D$, how can we find a $c$-confidence interval for the duration of the event?

In the comment thread in the link above, the following formula is proposed, where $g=1/\lambda$ is the mean interval for the Poisson process, and $Q^{-1}$ is the inverse of the regularized incomplete gamma function (InverseGammaRegularized in Mathematica):

$$\left\{g Q^{-1}\left(n,\frac{c+1}{2}\right),g Q^{-1}\left(n,\frac{1}{2}-\frac{c}{2}\right)\right\}$$

I suspect that the $n$ in the upper bound should be $n+1$, because without this, we get no upper bound if $n=0$.

Not understanding the math well enough to determine which formula is correct a priori, I wonder how we might go about testing a formula like this by simulation. My guess is that we should run many experiments, where in each experiment, we choose a random event $E$ and use a new Poisson process. If the formula is correct, the fraction of runs in which the event duration is within the bounds should be close to $c$. But how should we choose $E$?

Answer 1

Each Poisson arrival will see the system either in state $E$ or not-$E$ ($\bar{E}$). One might suspect that we can construct confidence intervals for the long run fraction of the time in state $E$ (LRF($E$)) by treating this as a sample from a binomial distribution for which we are estimating the probability $p$, but there are two problems with a straightforward approach. First, the constructed CIs only apply to the sample interval $D$. To see why this is a problem, imagine we sample over 10 seconds, 5 of which are in $E$ and 5 not. Increasing sampling frequency while holding $D$ fixed will cause our estimate to converge on 0.5 and our naive confidence intervals to shrink towards a width of 0, but the long run fraction of the time spent in $E$ might well not equal 0.5; $D$ is itself a sample from the long run.

Second, successive observations are not independent. Imagine the system alternates between one hour in $E$ and one hour in $\bar{E}$, forever. If the sample interval $D$ is, say, 10 seconds long with start time uniformly distributed over (0,2) hours, with high probability we will see exactly 1 sojourn in either $E$ or $\bar{E}$ and we will estimate LRF($E$) to be either 0 or 1. Our estimate will be very inaccurate, even with a sample size of 1000 over the 10 seconds. If, on the other hand, we have the duration of $D$ = 10 years with 1000 samples in $D$, we will see close to 500 sojourns in $E$ and 500 in $\bar{E}$ and our estimate will be close to 0.5. Our estimate, although based on the same sample size as the previous example's estimate, will be much more accurate. In either case, though, the estimator is unbiased, as its expected value is 0.5.

The other factor (besides sample size) that counts for constructing CIs for the LRF($E$) is evidently the number of distinct sojourns in $E$ and in $\bar{E}$ we see while in the steady state. The ideal case is when the sojourns are tagged, so our sample not only counts the frequencies of $E$ and $\bar{E}$ but the number of distinct sojourns into each state. Otherwise, if the mean duration of a sojourn in $E$ and $\bar{E}$ is much shorter than the mean inter-arrival time of our Poisson process, then we can assume almost all the observed $E$s are distinct sojourns, in which case the Binomial sampling approach gives only slightly too small confidence intervals. (This is pretty much the best case; it wastes the fewest Poisson samples.) If, on the other hand, the mean duration of a sojourn in $E$ and $\bar{E}$ is much longer than the mean inter-arrival time of our Poisson process, we might assume that a run of observed $E$s represents a single sojourn of the system in $E$, and likewise for $\bar{E}$.

Either way, by adding the sojourn counts for $E$ and $\bar{E}$ together, we get an estimated (or calculated) total number of sojourns $N$. Let $\hat{p}$ be the fraction of samples which saw state $E$; if $N\hat{p}$ and $N(1-\hat{p})$ are both large enough, say > 5 (rule of thumb), then we can construct an approximate CI using $\hat{p}$ and $\sqrt{\hat{p}(1-\hat{p})/N}$ as the mean and standard deviation in a Normal distribution, similar to what we would do with true Binomial sampling. Otherwise, we can construct an approximate CI using $N$ and, for the lower bound, $x = \lfloor N\hat{p} \rfloor$, for the upper bound, $x = \lceil N\hat{p} \rceil$, and pretending that we observed $x$ from a Binomial sample of size $N$.

Note that neither of these approximate CIs is likely to be any good when you're not in one of the two extreme cases described above.

Here's a little simulation that will illustrate the point. The duration of $E$ and $\bar{E}$ are both distributed Exponential(1), so the long run fraction of the time the system is in $E = 1/2$. I sample 1000 times at rates of every 0.01, 0.1, 1, and 10 time units, and repeat 1000 times, estimating $p$ (the long run average) each time. In the first case, it's easy to see we expect to see about 5 sojourns each in $E$ and $\bar{E}$ for an effective sample size of 10; in the latter, about 500 each, for an effective sample size of 1000. Here's the code and results:

MTBSamples <- c(0.01, 0.1, 1, 10)
EOn <- function(st) {min(which(E.Transition>st)) %% 2 == 1}
      # Odd = not in E, Even = in E

phat <- matrix(0, nrow=1000, ncol=length(MTBSamples))
for (i in 1:length(MTBSamples)) {
  for (j in 1:nrow(phat)) {
    E.Transition <- cumsum(rexp(20000,1))   #20000 = Lots more than needed
    SamplePoints <- cumsum(rexp(1000,1/MTBSamples[i]))
    phat[j,i] <- mean(sapply(SamplePoints, EOn))
  }
}

# Mean of estimates of long run occupancy fraction
colMeans(phat)
[1] 0.530586 0.502062 0.500477 0.499649

# Standard deviations of estimates of long run occupancy fraction
apply(phat,2,sd)
[1] 0.15686412 0.05153150 0.02290577 0.01681441

# (Estimated) effective sample size
0.25/apply(phat, 2, var)
[1]  10.15998  94.14438 476.48622 884.25297