Can Negative Binomial parameters be treated like Poisson?

Question

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

Answer 1

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!

Answer 2

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:

• The negative binomial occurs as a Poisson distribution compounded with a gamma distribution

  In this case the probability of success parameter, $p$ changes.

• The negative binomial occurs as a counting process where the interval/waiting time between events is geometric distributed.

  In this case the $r$ parameter changes.

---

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.

Answer 3

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