Number of events of type B given n events of type A in a Poisson process

I have a simple probability calculation I was working on that I came up with an answer for, but a question a colleague asked me led me to come up with a second approach - and a different answer. I put it aside at the time (a couple of years ago), but now it's niggling at me.

I think my original answer is correct because of other considerations, but I failed to identify a flaw in either approach; I'd appreciate responses that throw light on why one of them must be wrong.

Imagine two (independent) Poisson processes ('type A' and 'type B' say) happening along an interval; type A happens at rate $k_1 . \lambda$ and type B at rate $k_2 .\lambda$, where $k$'s are known (this could also be regarded as a partitioned Poisson process). However, we only see the type A events (we observe $n > 0$ events of type A). The idea is to try to get the distribution of the type B events given the observed type A. Putting the B events in lower case (to emphasize that we haven't seen them):

  | b  A  b b A   A   b   AA  b bb  A b |
0                                     1


Regarding type A as 'success' and B as 'failure' - recognizing that ignoring the actual times (which we don't actually observe) - and just considering the ordering of the events, it's immediately clear that the number of b's up to the $n^{th}$ A is NegBin$(n,p)$ where $p = k_2/(k_1 + k_2)$ (this being the 'count the failures' version of the Negative Binomial, rather than the 'count the trials' version).

The problem comes with how to deal with the number of events between the last A and the end of the interval.

One way to look at it is to now consider looking from the right end toward the left. The number of failures to the first success from the right is NegBin$(1,p)$, making the total number of type B events conditional on $n$ events of type A in the entire interval NegBin$(n+1,p)$.

The other way to look at it is to regard the whole interval as being on a circle. Now, between each consecutive pair of A's, the number of B's is NegBin$(1,p)$ and with $n$ such intervals the total number of B's is NegBin$(n,p)$.

(The first argument is the one the colleague's question led me to make)

By the memoryless property, it shouldn't matter where the interval starts; each time, the number of events to the next A should be NB$(1,p)$

Yet the answers cannot both be correct. Which is wrong?

• I agree with @Aniko : Since $N_2$ is independent of $N_1$ any knowledge about $N_1$ does not alter a probability about $N_2$. This is the answer to your question as it is stated. But you're actually interested in something else: as a statistician, you consider unknown parameters and you expect to predict something about $N_2$ by estimation. However, based on $N_1$ you can only estimate $k_1$, and you will not learn anything about $k_2$, hence you cannot predict $N_2$. Aug 5, 2014 at 14:57
• @StéphaneLaurent Both $k$s are known, not estimated. That's given in the question. Aug 5, 2014 at 21:14
• In your comments to @Aniko's answer you talked about the presence of unknown parameters. That's why I adressed this situation. Aug 5, 2014 at 22:05
• The unknown being the common $\lambda$. That is, $N_1$ contains information about $\lambda$ from which I want to form a prediction interval for $N_2$, which also depends on $\lambda$. Aug 5, 2014 at 22:47
• Ah indeed I missed $\lambda$. So indeed there's a negative binomial predictive distribution. I'm just waking up, I'll take a look later. Aug 6, 2014 at 6:58

Dare I say that both answers are wrong? Since $A$ and $B$ come from independent processes, the number of events of type $A$ in the $[0,1]$ interval has no effect on the number of events of type $B$ on the same interval. So the distribution of $B$s is still $Pois(k_2\lambda)$.

Here is a quick simulation study to demonstrate that:

R <- 10000
lbda <- 10
k1 <- 1
k2 <- 1.2
a <- rpois(R, lambda=lbda *k1)  # number of type A events
b <- rpois(R, lambda=lbda *k2)  # number of type B events
n <- 5
# select the samples where there are 'n' type A events
b.short <- b[a==n]
# plot results
plot(ecdf(b.short), main="Conditional distribution of B", xlim=c(0,max(b.short)))
curve(ppois(x, lambda=k2*lbda), col=2, add=TRUE)
curve(pnbinom(x, size=n, prob=k2/(k1+k2)), col=3, add=TRUE)
curve(pnbinom(x, size=n+1, prob=k2/(k1+k2)), col=4, add=TRUE)
legend(0, 1, c("Simulated","Poisson", "NB(n,p)", "NB(n+1,p)"),
col=1:4, lty=1, pch=c(19,NA,NA,NA), xjust=0, yjust=1)


Of course, this is not very helpful, because we do not know $\lambda$. We can, however, estimate it from the number of observed $A$s. From here, you would have to be specific about the estimate, and then one might be able to talk about the properties of $\hat{B}$.

EDIT: prediction interval

It seems that you are interested in the distribution of $B$ only for the purposes of a prediction interval. I think you can use your argument to build one. Let $B_n$ be the number of $B$ events before the $n$th $A$ event. Then $B_n \sim NB(n,p)$ and $B_A | A=n \sim NB(A,p)$, and for any given sample $B_n \leq B \leq B_{n+1}$. So it is not unreasonable to use the prediction intervals for either of them. At a first glance, the first should undercover, while the second overcover, but in practice discrete confidence intervals often overcover, so even the first one could be OK.

Just to follow up on the previous simulation:

coverage <- function(values, ints){
mean( (ints[,1] <= values) & (values <= ints[,2]))
}
pred1 <- cbind(qnbinom(0.025, size=a, prob=k1/(k1+k2)),
qnbinom(0.975, size=a, prob=k1/(k1+k2)))
pred2 <- cbind(qnbinom(0.025, size=a+1, prob=k1/(k1+k2)),
qnbinom(0.975, size=a+1, prob=k1/(k1+k2)))
# correct prediction interval
pred3 <- cbind(qpois(0.025, lambda=k2*lbda),
qpois(0.975, lambda=k2*lbda))

coverage(b, pred1)
[1] 0.9565
coverage(b, pred2)
[1] 0.9632
coverage(b, pred3)
[1] 0.9593


So pred1 undercovers with respect to the "correct" intervals (which overcover), and pred2 overcovers compared to it. Also note that a benefit of the $NB(n+1,p)$ interval is that it is defined even for $n=0$.

Just to summarize, neither of the arguments are correct, $B|A$ still has a Poisson distribution, but both of your distributions could be used for a reasonable prediction interval of $B$.

• Imagine both sets of observations were gaussian rather than Poisson. Your comment amounts to saying that a prediction interval for the next gaussian observation should use a normal distribution rather than a t-distribution. The presence of unknown parameters (lambda in the case of the poisson, sigma for the normal) alters the distribution. The uncertainty in lambda is what makes it negative binomial, the same way that uncertainty in sigma makes your prediction interval for a gaussian into a t-distribution. Jun 7, 2012 at 22:37
• Good analogy. But the distribution of the next observation, $Y_{n+1}$, is not $t$, it is normal! It is $\frac{Y_{n+1}-\bar{Y}}{\sqrt{S^2(1+1/n)}}$ that has a $t_{n-1}$ distribution, from which the prediction interval is derived. Note that it required estimates of the parameters, and with different estimates we would have obtained different intervals. Neither of your approaches mentions estimating $\lambda$, or estimating $B$ based on that, and the answer has to depend on that. Jun 8, 2012 at 13:27
• Take another look at the statistic you mentioned that has a $t$ distribution. It only involves an estimate of $\mu$ in the sense that you can say 'that second term on the numerator estimates $\mu$' - but the expression itself doesn't - $\mu$ has been 'removed' by subtracting the two terms involving it (which is why you can work out the distribution without reference to it). In the Poisson, the quantities in the estimate of $\lambda$ are in the Negative Binomial expression. A different estimate of $\lambda$ corresponds to a different interval. I can make $\lambda$-hat explicit if it helps. Jun 8, 2012 at 23:08
• So, for example, the 'wrap it round the circle' version is NegBin$(k_1 \hat{\lambda}, \frac{k_2}{k_1 + k_2 })$ which has expectation $k_2 \hat{\lambda}$ and variance that is about what you'd expect - the sum of the variance of the Poisson at $\lambda = \hat{\lambda}$ plus the variance of the estimate of the mean. Jun 9, 2012 at 5:20
• I added some further discussion to my answer. By the way, I love your question, it's a wonderful brain-twister. Jun 11, 2012 at 13:40

Sorry for having taken so much time to correctly read the statement of the problem. The Poisson processes just hide the problem. Basically you have two independent Poisson variates $X_1 \sim {\cal P}(\lambda k_1)$ and $X_2 \sim {\cal P}(\lambda k_2)$ with known $k_1, k_2$ and unknown $\lambda$, and based on a realization $x_1$ of $X_1$ the goal is to predict the realization $x_2$ of $X_2$.

A Bayesian approach with the Jeffreys/Bernardo prior provides a prediction interval enjoying very good frequentist properties. Assume the Jeffreys prior $\Gamma(\frac12,0)$ on $\lambda$. Then after observing $x_1$ the posterior distribution on $\lambda$ is $\Gamma(x_1+\frac12, k_1)$. Integrating the distribution of $X_2$ with respect to $\lambda$ over this posterior distribution yields the ${\cal P}{\cal G}(x_1+\frac12, k_1/k_2)$ (Poisson-Gamma) distribution as the posterior predictive distribution of $x_2$, which is the same as the negative binomial distribution ${\cal N}{\cal B}(x_1+\frac12, k_1/(k_1+k_2))$. Taking quantiles of this distribution provides a prediction interval about $x_2$.

• Thanks, yes, a Bayesian approach is easy for this problem (I've done it on this problem enough times I can do that one in my sleep). However, I'm trying to arrive at a frequentist prediction interval. The Poisson process is a way to try to obtain what corresponds to a pivotal quantity. I assume that if one takes a flat prior one should get the same result as the frequentist prediction interval (as that's often the case), but I don't see how to prove that for this problem. Aug 6, 2014 at 9:49
• @Glen_b There's a usual method for a frequentist prediction interval for a Poisson sample, but I don't remember it (it is in a well-known book about statistical intervals). There should be no difficulty to apply it here. $k_1$ and $k_2$ are like sample sizes here. Imagine $k_1$ and $k_2$ are integers, then the problem is equivalent to the situation of a ${\cal P}(\lambda)$ sample of size $k_1$ and you want to predict the sum of a ${\cal P}(\lambda)$ sample of size $k_2$. Aug 6, 2014 at 9:54
• Thanks Stéphane. I had presumed it would be an already solved problem, but I hadn't seen it anywhere. (By the way, it's not that I have anything against a Bayesian approach for the problem this is a simplified version of -- it's ideal as far as I'm concerned -- but the choice is not mine to make.) Aug 6, 2014 at 9:55
• Frequentist method in section 3 Aug 6, 2014 at 10:24
• Yes, thanks again. I've actually read the Bain and Patel paper before (which is in the references of that paper), which covers backing endpoints out form a binomial (as in Sec 3.1). Unfortunately, that's not of direct help in my problem; I was trying to compute the distribution of a pivotal or approximately pivotal quantity in order to compute convolutions of several of these -- in order to predict a sum of independent terms of this form. (Given the post was already very long, I felt adding in all that details would make it overly daunting.) ... (ctd) Aug 6, 2014 at 10:58