I have a problem where I need to compute the likelihood of a given sequence of events (in time) given a Poisson process with a particular lambda (events per second). I only need it to within a constant, for comparison to other possible sequences. A likelihood ratio, for example:

p(sequence_1 | lambda) / p(sequence_2 | lambda)

These probabilities could be computed by computing the value of the Poisson distribution over the full time span. For example, if the whole interval is 10 seconds, I could compute the value of the Poisson distribution with 10*lambda and the total number of events.

However, I could also break the full interval down into subsets, such as 1 second intervals. I think the total probability for a sequence ought to be related to the product individual Poisson probabilities for each sub-interval, but I think I must be missing a normalizing factor. Furthermore, at an intuitive level, if all the events were bunched up in time, analyzing over the full time may appear likely while looking at the sub intervals is unlikely.

Is there a way to do this where the result is not a function of the intervals used for processing?

EDIT: For clarification, my goal is to check a sequence against the Poisson process and see how well it matches. I need the score to be in units of probability such that if the score is twice as high, it is twice as likely to be from the specified Poisson.

  • 1
    $\begingroup$ What you're calling a likelihood ratio is actually a probability ratio: you are using the pipe (conditional) operator, but this isn't a Bayesian treatment. I'm sure most agree a likelihood ratio would be written as $L(\lambda_1, X) / L(\lambda_2, X)$ with $X$ the same observed data. The LR is used to compare parametrizations for observed data, not to combine or parcel out values. $\endgroup$
    – AdamO
    Dec 16, 2021 at 23:34
  • 1
    $\begingroup$ The occurrence of "$N$ events in an interval of length $L$" (having a Poisson probability) is not the same as the conjunction of the events "$N_1$ events in the first interval of length $1,$" "$N_2$ events in the second interval of length $1,$" ..., "$N_L$ events in the last interval of length $1,$" where $N_1+N_2+\cdots N_L=N.$ The latter is far more detailed and therefore has a different (much smaller) likelihood. When you account for all the (multinomial) combinations of how the $N$ events can fall into the $L$ unit intervals, you get precisely the original Poisson probability. $\endgroup$
    – whuber
    Dec 17, 2021 at 0:24
  • $\begingroup$ @whuber, that all makes sense. But in the end, is there a "right" or best way to see the goodness of fit of a sequence to a particular Poisson (lambda)? $\endgroup$
    – Mastiff
    Dec 17, 2021 at 15:17
  • 1
    $\begingroup$ That's a very helpful reformulation of your question. Because most people will overlook it (I did until just now), please include it in an edit. There definitely are good solutions, but the simplest (and arguably best) ones don't rely on computing likelihoods of individual sequences. Assuming you have sufficient data (a practical minimum would be ten events), a chi-squared test based on a predetermined number of subdivisions of the interval will work well, because it conforms to all the requirements laid out in my example at stats.stackexchange.com/a/17148/919. $\endgroup$
    – whuber
    Dec 17, 2021 at 16:15
  • $\begingroup$ I tried to keep things simple in the question, but this often backfires since people innocently work around the actual problem. The application here is a hypothesis test of possible sequences. The "Poisson-ness" of the sequence is one component among many I am scoring, in units of probability. My goal is to be able to score them all such that their scores are proportional to their relative probabilities. $\endgroup$
    – Mastiff
    Dec 17, 2021 at 17:26

2 Answers 2


Consider a chi-squared test.

Before examining the data, you may partition the time interval (which we can, by establishing a suitable origin and unit of time, assume to be $(0,1]$) into $k\ge 2$ subintervals $\mathcal{I}_j=(t_{j-1}, t_{j}]$ with $0=t_0 \lt t_1 \lt \cdots \lt t_k = 1$ and count the numbers of events in each subinterval.

Under the null hypothesis of a homogeneous Poisson process of intensity $\lambda,$ the expected number of events in $\mathcal I_j$ is $\lambda(t_{j} - t_{j-1}).$ The maximum likelihood estimator of $\lambda$ is $n$ itself. The alternative hypothesis of interest is any departure from this condition of homogeneity. Thus, this situation conforms exactly to the description of the chi-squared test at https://stats.stackexchange.com/a/17148/919. The reference chi-squared distribution for computing the p-value has $n-1$ degrees of freedom because one parameter has been estimated.

Given an array x of the observed times in $(0,1]$ of any continuous point process, here is an R function to carry out the chi-squared test.

PPGOF <- function(x, k, ...) {
  n <- length(x)
  if (missing(k)) k <- ceiling(sqrt(n))
  if(isTRUE(n < 5*k)) k <- ceiling(n/5)
  if(isTRUE(k < 2)) k <- 2
  chisq.test(tabulate(ceiling(x*k), k), ...)

The variable k is the number of equal-width subintervals to create. (I leave coding the more general test, with variable length intervals, to anyone who needs this: the code is fussier but not essentially different.) This solution merely tabulates the counts within each subinterval and passes that summary on to the built-in chisq.test function. The result is an object that includes the chi-squared test statistic and its p-value, which can be extracted via $p.value (as shown in the code at the end of this post).

How effective is this? Here is an example with $n=100$ events. The top row displays (a) the (constant) intensity function, plotting the relative rate $\lambda(x)$ against the time $x;$ (b) one realization of this process; and (c) the sampling distribution of chi-squared p-values (based on ten thousand independent realizations). Because this is the null hypothesis, ideally that distribution is uniform. Clearly it approximates uniformity closely.


The bottom row is a similar array of displays for a non-constant intensity function: that is, a non-homogeneous Poisson process. To my eye, at least, there is no qualitative difference between the illustrated realization of the process and the preceding realization of a homogeneous process: evidently, the difference is subtle. Nevertheless, as the bottom right plot shows, the sampling distribution of p-values is markedly different: they tend to be higher than before, especially in the low range. In other words, this test has appreciable power to determine that this process is not homogeneous.


This reference chi-squared distribution will work well for large $k$ and/or large $n.$ For small $k,$ to rely on the chi-squared distribution for computing p-values, you need to respect the standard rule of thumb that all of the expected counts are $5$ or greater. If this is not possible, you can resort (as is usual) to permutation testing. For small test sizes and smallish $n$ you will discover that the discreteness of the sampling distribution of the statistic might make the test more rigorous than desired, so make sure to study this possibility before proceeding in such cases.

If you anticipate certain forms of inhomegeneity might occur, it is best to partition time into intervals with approximately equal expected counts. The flexibility afforded by choosing the numbers and widths of these intervals makes the chi-squared test especially appropriate and powerful in this application. However, one thing you do not want to do is experiment with choices of $k$ and the $t_j:$ that form of "data snooping" would yield unrealistically small p-values. It would be OK, though, to use a preliminary dataset (perhaps a random subset of the data) to choose $k$ and the $t_j$ and then apply those to the rest of the data.

Code example

The following code creates the figures and illustrates how PPGOF can be used.

# Inhomogeneous Poisson process.
# f: [0,1] --> R will be reduced modulo 1.
rPP <- function(n, f=identity) {
  x <- cumsum(rexp(n+1))          # Use exponential waiting times
  f(x[-(n+1)] / x[n+1]) %% 1      # Apply `f` to these homogeneous Poisson events

# Describe the simulation study.
n <- 100     # Sample size to simulate
n.sim <- 1e4 # Number of iterations
hypotheses <- list(`0` = identity,                          # Null: homogeneous process
                   `A` = function(x) x * (25 + sin(4*pi*x)))# Alternative process

# Carry out the study.
mai <- par("mai"); par(mfrow=c(length(hypotheses), 3), mai=c(.5, .4, .5, .1))
for (h in names(hypotheses)) {
  # Compute and plot the intensity. This works for increasing functions f:[0,1] --> R+.
  f <- function(x) {g <- hypotheses[[h]]; g(x) / g(1)}
  x <- seq(0, 1, length=201)
  y <- diff(x) / diff(f(x))  # Numeric derivatives (intensities)
  plot(x[-1], y, type="l", ylim=c(0, max(y)), xlim=0:1, lwd=2, col="Blue",
       xlab="x", ylab="Intensity", main="Intensity Function")
  abline(h = 1, lty=3) # The reference line is the homogeneous intensity
  # Plot an instance of this process.
  PP.plot(rPP(n, f), col=gray(0, alpha=.25), main=bquote(paste("Instance of ", H[.(h)])))
  # Estimate the sampling distribution and plot it.
  H <- replicate(n.sim, PPGOF(rPP(n, f))$p.value)
  plot(ecdf(H), xlim=0:1, xlab="p value", lwd=2,
       main=bquote(paste("Sampling Distribution under ", H[.(h)])))
  abline(0:1, col=hsv(0, alpha=0.5), lwd=2) # The reference line is a uniform distribution
par(mfrow=c(1,1), mai=mai) # Restore the original plotting environment

This is a basic probability result: if $Y_1, Y_2, \ldots, Y_n \sim \text{Poisson}(\lambda)$ then $\sum_{i=1}^n Y_i \sim \text{Poisson} (n\lambda)$. So differing from the premise of your question: the "result" is never a function of the intervals used for processing.

The combined and the separate analyses produce the same "result" whether a prediction, inference, estimation, etc. Dividing up a Poisson sample into the independent contributions makes no impact on a final evaluation of probability or likelihood like we would intuitively expect.

  • $\begingroup$ In practice it does, or perhaps I'm missing your point. If I compute the probability of the total events over an hour and compare it to the product of the probabilities of the number of events in each minute, the answers will be different. A conceptually similar idea with binomial: the odds of getting 1 million heads out of two million flips is relatively high, but if all 1 million heads came in the first million flips and I scored the first and second halves separately, the probability would be infinitesimal. $\endgroup$
    – Mastiff
    Dec 16, 2021 at 23:48
  • $\begingroup$ Should you be looking at the sum of the probabilities of the number of events in each minute over an hour? $\endgroup$ Dec 17, 2021 at 0:10
  • 1
    $\begingroup$ @Geoffrey Such a sum could easily exceed $1$ when the rate is approximately one per minute or the hour is prolonged. $\endgroup$
    – whuber
    Dec 17, 2021 at 0:17
  • $\begingroup$ "One million heads out of two million flips" comprises an astronomical number of sequences, whereas "all million heads came in the first million flips" is a unique sequence. That's why the former has a hugely greater chance than the latter. $\endgroup$
    – whuber
    Dec 17, 2021 at 0:27
  • 1
    $\begingroup$ A chi-squared test will address that question nicely. $\endgroup$
    – whuber
    Dec 17, 2021 at 17:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.