Estimating event duration by sampling with a Poisson process

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$?

• To help you clarify your thinking, please consider this: why do you consider $E$ to be "random"? In your considerations, also ponder whether the set $E$ itself matters in any of the calculations or simulation, or whether all that matters is the measure of $E$.
– whuber
Jun 24, 2012 at 14:30
• Having read the link, it's not clear to me whether you really mean "duration" or, alternatively, "long-run fraction of the time the system is in $E$". The two are different; if the latter, there's a much simpler way of calculating a confidence interval from a sample of Poisson arrivals. Jun 24, 2012 at 15:09
• whuber, I suspected that all that matters is the measure, and so ran experiments where, for each sampling time, the result was a Bernoulli trial with success probability equal to the measure of E divided by the duration of the experiment. I feel that E must be chosen "randomly" each time due to my understanding of the definition of a c-confidence interval. So for example, using Bernoulli trials, I thought it necessary to choose the measure of E randomly. Not sure what you are getting at.
– Tom
Jun 24, 2012 at 18:44
• jbowman, I think we want "long-run fraction of the time the system is in $E$", like you say. I see the formula above as providing an estimate for this given the currently available data (after dividing by $D$). What's the difference? In any case getting estimates of total event duration during an experiment is also useful, but long-run fraction is more useful.
– Tom
Jun 24, 2012 at 18:48
• Tom, the randomness is in the Poisson process. It matters not whether E is fixed or random, provided you do not change its measure. For a simulation study, then, you would create a sequence of E's of different measures and resample the Poisson process multiple times.
– whuber
Jun 25, 2012 at 14:18

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

• I had a hard time getting through the first paragraph: it looks like something got garbled during the posting. Do you think you could clean that up?
– whuber
Jun 25, 2012 at 14:19
• @whuber ... or I was channeling James Joyce when I wrote it. I'll redo the post later today (it's midday my time.) BTW +1 for the tact! Jun 25, 2012 at 18:44
• Thanks for the detailed answer. This is much more complicated than I expected, and still (at least!) a little above my head.
– Tom
Sep 20, 2012 at 4:22
• I'm left wondering whether the model you use for the events is good for the specific application — for example, after thinking about it more, it seems hard to define the "long-run fraction of the time" in this case, since the events are all finite (i.e. the system is in state "dead" almost all of the time). I suppose you have to model the data as if your current behavior were extended forward through all time? Going to go ahead and mark this as the best answer now, though.
– Tom
Sep 20, 2012 at 4:28