8
$\begingroup$

I have a count process that I'd like to model with a Poisson process. Data is measured every 30 minutes, and with a poisson distribution I can easily measure the probability of a given count of events being anomalous in different time periods using a fitted value of lambda, i.e. "is the number of events we've seen in the 30 min anomalous? What about the last hour? Is the number of events we've seen in the last 1.5 hours anomalous?", etc.

The problem is that my data is overdispersed, and definitely is described well by a negative binomial distribution. I'm choosing to use the parameters $(\mu, \alpha)$ since that's what PyMC3 uses, where $\mu$ is equivalent to lambda from the poisson distribution.

Is there a way to utilize the negative binomial parameters in the same way as the poisson rate parameter where I see if an event count is anomalous in some time period t (where I can extend t to different time periods)?

$\endgroup$
1
  • $\begingroup$ Is there a way to utilize the negative binomial parameters in the same way as the poisson rate parameter where I see if an event count is anomalous in some time period t (where I can extend t to different time periods)? In what way did you do the anomaly detection with Poisson distributed data? And what do you mean by 'the same way'? Do you mean to use the same principle of anomaly detection or do you want to use the exact same method, including the use of a Poisson distribution as assumption to do fitting (only lambda and not $\mu$ and $\alpha$? $\endgroup$ – Sextus Empiricus Sep 19 '20 at 18:14
7
+50
$\begingroup$

I put up some code to perform this task in PyMC3, since you mentioned it in the question. The first part, which you seem to already be familiar with, would be fitting the model to get a posterior distribution on the parameters:

import pymc3 as pm
import numpy as np

# generating simulated data data for a week
data = pm.NegativeBinomial.dist(mu=3, alpha=1).random(size=7*24*2)

# defining the model and sampling (MCMC)
with pm.Model() as model:
    alpha = pm.Exponential("alpha", 2.0)
    mean = pm.Exponential("mean", 0.2)
    obs_data = pm.NegativeBinomial("obs_data", mu=mean, alpha=alpha, observed=data)
    trace = pm.sample()

# plotting the posterior
pm.traceplot(trace)
pm.plot_posterior(trace)

Now we get to the part on which you seem to be struggling. We can use this nice property: when two random variables, $X$ and $Y$ have negative binomial distributions with the same overdispersion parameter, then $X+Y$ also has negative binomial distribution, with mean $\mathbb E[X]+\mathbb E[Y]$ and the same overdispersion parameter as $X$ and $Y$. You can find proofs for this property here.

Assuming that the negative binomial parameters are fixed (formally, assuming your stochastic process is in the class of Lévy processes, in which Poisson processes are included), that implies that if you want to know the distribution for the number of events in a whole hour or a whole day, you just have to adjust the mean, like you would do with a Poisson process.

For example, to find out how atypical it would be to find more than 200 events in a single day, we could use the following:

np.mean(pm.NegativeBinomial.dist(mu=48*trace["mean"], alpha=trace["alpha"]).random(10**4)>200)

Let's break this line of code down for a bit. When we use pm.NegativeBinomial.dist(mu=..., alpha=...), we are invoking the PyMC3 implementation of the negative binomial with a specific set of parameters (we could use the Numpy implementation as well, but they are parameterized differently, so it is less error prone to stick to PyMC3).

We then use the parameters we sampled from the posterior: alpha=trace["alpha"] for the overdispersion, and mu=48*trace["mean"] for the mean (we multiply by 48 to adjust this mean to reflect 24 hours instead of half an hour).

Finally, we sample many instances from this distribution and compare them to the value we are interested in (.random(10**4)>200), then finding the probability of new samples from our negative binomial process exceeding it (by applying np.mean to the resulting array of booleans). The result is the probability that your model would generate a day with 200 events or more.

A few caveats here:

  • if your model allows overdispersion to change with time, none of this will work
  • if your model allows the Poisson rate to change with time as some function $\lambda(t)$, this will have to be somewhat adapted. Instead of multiplying the rate by some number, you will have to integrate $\lambda(t)$, making this a bit more complicated.

EDIT: I am editing to address the comment by @J Does asking about day of week effects. So, let us first generate some data with strong day of week effects:

# how many weeks of data are available?
WEEKS = 5
# how many observations are available per day?
OBS_PER_DAY = 24*2

data = pm.NegativeBinomial.dist(mu=[2,3,1,2,5,9,7]*5, alpha=1).random(size=OBS_PER_DAY).T.flatten()

Now, one way we can get around that is by having 7 different means, instead of a single one. The PyMC3 model can be written as:

with pm.Model() as model:
    alpha = pm.Exponential("alpha", 2.0)
    mean = pm.Exponential("mean", 0.2, shape=7)
    day = np.arange(WEEKS*7*OBS_PER_DAY)//OBS_PER_DAY%7
    obs_data = pm.NegativeBinomial("obs_data", mu=mean[day], alpha=alpha,
        observed=data)
    trace = pm.sample()

The variable day here associates each observation to the day of the week it came from. Now, we have a model that allows for day of week effects. How can we check if having more than 500 events on a Friday is atypical? The procedure is similar to the homogeneous case:

friday = 4 # assuming the week starts on monday
np.mean(pm.NegativeBinomial.dist(mu=48*trace["mean"][:,friday], alpha=trace["alpha"]).random(10**4)>500)

OK, now what if we want to check if 3000 events in a week is an atypical event? The expected count of events for a week is 48*sum(mean), so we do this:

np.mean(pm.NegativeBinomial.dist(mu=48*trace["mean"].sum(axis=1), alpha=trace["alpha"]).random(10**4)>3000)

Notice that we did not need any fancy integration, since this day-of-week effect makes $\lambda(t)$ a piecewise constant function. (hooray!). You will need no integrate the Poisson rate when its functional form is a bit more complicated: for instance, if $\lambda(t)$ is a polynomial, an exponential, a function sampled from a gaussian process, etc. Unfortunately, it seems to be hard to find resources on this specific topic on the Web... Perhaps I will add something addressing this issue to this answer when I find the time.

Hope I was helpful!

$\endgroup$
2
  • 1
    $\begingroup$ This is exactly what I was looking for thank you, appreciate how you included the proofs as well. Question about that second caveat - If I were to have a negative binomial process like this but the poisson rate had day of week seasonality, do you have any resources that might help me to learn that integration process? It seems pretty useful to know so might as well! $\endgroup$ – J Doe Sep 19 '20 at 20:54
  • $\begingroup$ @JDoe I edited the answer to address day of week effects. I couldn't find much resources on the integration of Poisson rates $\lambda(t)$ for this purpose, also no implementations in PyMC3... I will try to add something about it in this answer eventually. But the particular case of day of week effects is easier and is now addressed there $\endgroup$ – PedroSebe Sep 20 '20 at 2:44
3
$\begingroup$

The negative binomial can be treated like Poisson, but it is ambiguous how to treat it. It will depend on the underlying process that causes the overdispersion. This can occur in different ways.

Below I will describe two ways:


1. Compound distribution

You can view the negative binomial distribution as a Poisson distribution compounded with a gamma distribution.

If

$$Y \sim Poisson(\lambda=X)$$ where $$X \sim Gamma(\alpha,\beta)$$

Then $$Y \sim NB(r=\alpha, p = (\beta+1)^{-1})$$

With a Poisson process, if you consider a larger time interval, then the distribution of the number of events relates to a Poisson distributed variable with a larger rate coefficient.

For instance, the Poisson rate in the compound distribution is scaled with a factor $c$.

$$Y_c \sim Poisson(\lambda=cX)$$

This is similar to scaling the rate of the gamma distribution.

$$cX \sim Gamma(\alpha,\beta/c)$$

So the compound distribution becomes

$$Y_c \sim NB(r=\alpha, p = (\beta/c+1)^{-1})$$


2. Counting process

You can view the negative binomial distribution as occurring in a counting process where the waiting time between events is i.i.d geometric distributed.

If you consider the ordered sequence of events $1,2,...,k,k+1,...$ where the time between events follows a geometric distribution:

$$t_k-t_{k-1} \sim Geom(p)$$

Then the number of events within an interval of length $t$ follows a negative binomial distribution with $r=\lfloor t \rfloor$ and $p=p$

$$N_{\text{events within $t$}} \sim NB(\lfloor t \rfloor, p)$$

In that case, the increase of the time period $t$ over which the counting process is performed corresponds to an increase of the parameter $r$ in the negative binomial distribution.

This case corresponds to the answer by PedroSebe.


So it will depend on what sort of process you have that generates the negative binomial distribution of counts.

$\endgroup$
4
  • $\begingroup$ Thanks for this! I'm a little confused so I'm going to try and state it back to you. So say I fit NB(r_x, p_x) to a time series of counts (measuring frequency is every x=30 minutes), and then I have new data come in, its a count of events in the last y=45 minutes. I can determine the probability of that many events occurring in 45 minutes given the data is from a NB(r_x, p_x) distribution by calculating the PMF of NB(r_x, p_y), where p_y = 1 / ( (1/p_x - 1) * 30/45 + 1) )? And I don't need to scale the parameter r_x to adjust to time y=45 minutes? $\endgroup$ – J Doe Sep 19 '20 at 20:22
  • $\begingroup$ @J Doe , it will depend on the origin of the overdispersion. I guess that PedroSebe's answer is more typically found when you have a situation that the source of overdispersion is a waiting time that is not exponential distributed (although I am wondering what sort waiting time distribution corresponds to a negative binomial distribution for number of events and imagine that it is only approximately negative binomial distributed). The situation I used, where r parameter is constant, is when you consider a waiting time that is exponential distributed but the rate varies over a long timescale. $\endgroup$ – Sextus Empiricus Sep 19 '20 at 23:01
  • $\begingroup$ @SextusEmpiricus using a parameter r != 1 means the waiting time distribution is a Gamma instead of an Exponential. $\endgroup$ – PedroSebe Sep 20 '20 at 2:49
  • $\begingroup$ @PedroSebe I mean, for number of events distributed by the Poisson distribution the corresponding waiting time for time between events is exponential distributed (or gamma distributed when you wait for the n-th event). For the binomial distribution you could see the waiting time as a geometric distributed variable. But what would be the corresponding waiting time distribution, if there is one, for a negative binomial distribution? $\endgroup$ – Sextus Empiricus Sep 20 '20 at 6:42
1
$\begingroup$

This is how I would do it in R. If correct, it should be easy to translate to python.

First estimate the parameters of the negative binomial distribution that best fit a given training dataset. Then map the new data to the distribution function with those parameters.

library(MASS)

set.seed(1234)
data_stream <- rnbinom(n= 1000, size= 1, mu= 10)

params <- fitdistr(x= data_stream, densfun= 'negative binomial', lower= c(1e-9, 0))
params
          size           mu     
   0.96289937   10.02900002 
 ( 0.04719405) ( 0.33835666)

new_time_point <- 30

pnbinom(new_time_point, size= params$estimate[1], mu= params$estimate[2])
0.94562 # This is how extreme the new data is
$\endgroup$
1
  • $\begingroup$ Thanks for the reply! I already know how to fit negative binomial parameters to data and how to make point estimates using distribution functions. Let's say each observation in data_stream was the number of counts seen in a given 30 minute interval. And let's say I've only seen 1 count in the past 80 minutes. I'd like to know how extreme is it to only have seen 1 count in the past 80 minutes? The difference between what I'm asking for and your answer is extending point estimates (with pnbinom for instance) to different time periods than the measuring interval $\endgroup$ – J Doe Sep 14 '20 at 15:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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