Confidence interval over multiple Poisson distributions

Question

There are five different periods (one following each other) of unequal lengths where one specific event is happening such that it follows Poisson distributions (e.g. meteorite falls). It is always the same type of event, just the expected average time between event is changing from time to time. In reality there are many more such periods.

The data look like this:

• in the 1st period (100 hours) the average interval for event is 500 hours
• in the 2nd period ( 80 hours) the average interval for event is 700 hours
• in the 3rd period (130 hours) the average interval for event is 300 hours
• in the 4th period ( 50 hours) the average interval for event is 800 hours
• in the 5th period (400 hours) the average interval for event is 900 hours

Lambda:

• 1st period lambda: (1/500)*100 = 0.200
• 2nd period lambda: (1/700)*80 = 0.114
• 3rd period lambda: (1/300)*130 = 0.433
• 4th period lambda: (1/800)*50 = 0.062
• 5th period lambda: (1/900)*400 = 0.444

Is there a way to estimate how many events happened and what is the confidence interval for the estimate?

I understood from How to calculate a confidence level for a Poisson distribution? and here: http://ms.mcmaster.ca/peter/s743/poissonalpha.html that I can use (x is number of events):

• ( qchisq(0.025, 2*x)/2, qchisq(0.975, 2*(x+1))/2 )

or maybe some other method that will be more suitable for low lambda.

I don't understand how I can add all the five distributions together to say how many events were expected.

Answer 1

Because (a) the number of events in any $T$ hour period with a rate of one event every $H$ hours follows a Poisson distribution with parameter $T/H$ and (b) the number of independent Poisson random variables of rates $\mu_1,\mu_2,\ldots,\mu_n$ in non-overlapping intervals is a Poisson random variable of rate $\mu=\mu_1+\mu_2+\cdots+\mu_n,$ any question about the probability distribution of the number of events is readily answered. Let's write $F$ for this distribution function and (therefore) $F^{-1}$ for its quantile function.

For instance,

• The expected number of events is $\mu.$

• The variance of the expectation is $\mu$ and so its standard deviation is $\sqrt{\mu}.$

• For any tolerance level $q$ between $0$ and $100%,$ choose nonnegative values $\alpha_l$ and $\alpha_u$ totaling $1-q.$ Then with at least probability $q$ the number of events will be between $F_{l} =F^{-1}(\alpha_l)$ and $F_{u}=F^{-1}(1-\alpha_u).$ The actual probability is $F(F_{u}) - F(F_{l}-1).$

Often one chooses a "symmetrical" interval in the sense that $\alpha_l=\alpha_u = (1-q)/2.$

To see how simple this is, let's use the data in the question as an example.

1. The time intervals are $100,80,130,50,400$ minutes and the rates are the reciprocals of $500,700,300,800,900$ hours, whence $$\mu = \frac{100}{500} + \frac{80}{700} + \cdots + \frac{400}{900} = 1.2545\ldots\,.$$

2. The standard deviation is $\sqrt{1.2545\ldots} \approx 1.12.$ You can use this immediately (using only mental arithmetic) to bound the chances of large numbers of events. For instance, Cantelli's Inequality implies the chance of $5$ or more events, which is $(5-1.2545)/1.12 \approx 3.3$ standard deviations beyond the mean, cannot exceed $1/(1+3.3^2) \approx 1/12 \approx 8\%.$

3. If, say, you wanted a $q=95\%$ symmetric interval in which the number of events is likely to fall, set $\alpha_l=\alpha_u=2.5\%$ and locate the $0.025$ and $0.975$ quantiles of the Poisson$(1.2545)$ distribution at $0$ and $4,$ respectively. The actual chance that the number of events will be between $0$ and $4$ inclusive is $99.1\%.$

---

To perform the calculations in the example I used the following commands in the statistical calculator R.

#
# Specify the problem.
#
intervals <- c(100, 80, 130, 50, 400)
averages <- c(500, 700, 300, 800, 900)
q <- 0.95
#
# Compute some properties of the distribution of the total number of events.
#
mu <- sum(intervals / averages)                    # Expected number
alpha.u <- alpha.l <- (1 - q)/2                    # Upper and lower tolerances
tl <- qpois(c(Lower=alpha.l, Upper=1-alpha.u), mu) # Limits
coverage <- diff(ppois(tl + c(-1,0), mu))          # Probability within limits 
names(coverage) <- "Coverage"
#
# Display the answers.
#
(tl) # The limits
(signif(c(Expected=mu, SD=sqrt(mu), `Nominal coverage`=q, coverage), 3))