# Establishing the Probability Distribution that Governs a Random Process

I have a experiment in which a random event that can occur $N$ times in some set period of time $T$ (where $N \in$ $\mathbb{N_0}$ but is typically less that 10 and never exceeds 100). Here I will refer to a single experiment as a "session", and for a single session we know that the probability of $n$ events occurring in some remaining time $T_{R} = T - t_{i}$ is approximately Poisson distributed, but that this distribution is not "quite right".

In order to attempt to model the exact way these events are distributed during a "session", I have run ~60,000 sessions, each one with $T = 60$ minutes. From these "sessions" I have gathered the data of how many "event" occurred in the time T. This looks like

We can model the the probability of a number of events occurring in the period $0 \leq t < T$ using a Poisson distribution (which is providing the probability of $n$ events occurring in the time-remaining $t$).

To model the events we let $t$ be the time remaining and $X_t$ be the total number of events to occur in the remaining time $t$. The model implies the estimates $P(X_t \leq c)$. So under our assumption $X_t\sim \operatorname{Poisson}(\lambda_{t})$ we have: $$P(X_t \leq c) = e^{-\lambda}\sum_{k=0}^c\frac{\lambda_t^k}{k!}\,.$$

This works fairly well during simulation (where we have an implied $\lambda$ from another data source), but we have a belief that the underlying process can be more accurately modeled using a Negative Binomial distribution rather than Poisson. My problem is that I am not a statistician and I don't really know how I can go about extracting the relevant parameters for the Negative Binomial using the data above. I have access to R and am happy to write the necessary code to extract these parameters, I am just unsure about what this process/method of extraction is, can anyone outline this for me?

Better still, if we switch to use Negative Binomial distribution so that we assume now that the probability of a number of events occurring in the period $0 \leq t < T$ is $X_t\sim NB(r, p)$ and we have $$P(X_{t} \leq c) = p^{r}\sum_{k = 0}^c (1 - p)^{k}\binom{k + r -1}{r - 1},$$ can anyone provide an R script that can be used to extract my $r$ and $p$?

The data above is:

╔═══════════════╦══════════╗
║ Session Count ║ Events   ║
╠═══════════════╬══════════╣
║  26069        ║ 0        ║
║  30175        ║ 1        ║
║  18997        ║ 2        ║
║  8136         ║ 3        ║
║  2934         ║ 4        ║
║  820          ║ 5        ║
║  250          ║ 6        ║
║  54           ║ 7        ║
║  16           ║ 8        ║
║  4            ║ 9        ║
║  7            ║ 10       ║
║  1            ║ 11       ║
║  0            ║ 12       ║
║  1            ║ 13       ║
╚═══════════════╩══════════╝

• Keep in mind that the Poisson distribution is equivalent to the NB distribution in the limit that $r \rightarrow \infty$ and $p \rightarrow 0$ such that $r p = \lambda$ is a constant. This has a nice interpretation when your "samples" are over a continuum of time instants. A "sample" can be interpreted as an instant of time, with infinitesimal probability of success, but there are infinite samples. I can't think of a good interpretation of $r$ and $p$ in this case, if they're finite. I don't know if interpretability is what you're after, but in case it is, I thought this would be nice to know. – Bridgeburners Feb 8 '18 at 21:44
• Isn't this answered at stats.stackexchange.com/questions/164934/… and stats.stackexchange.com/questions/80026/…? Some useful commentary also appears at stats.stackexchange.com/questions/173450/…. If any of these answers your question, let us know before the bounty runs out and we can return it to you. – whuber Feb 8 '18 at 22:05
• @whuber thanks for your suggestions, you time is most appreciated. I will check later today and get back to you. – MoonKnight Feb 9 '18 at 9:18
• @MoonKnight Is simultaneous estimation of $r$ and $p$ in this wiki perhaps what you are looking for? – Jim Feb 11 '18 at 20:38
• I know this gets to extremely unlikely events, but in both a Poisson and NB process technically there is a probability that the number of events observed in a period is any value up to infinite... so if you are interested in imposing a hard 100 count ceiling, you might want to consider a truncated distribution of some form. – Alexis Feb 13 '18 at 15:52

Given some data and the general shape of a distribution, it's possible to fit the parameters of that distribution using the principle of Maximum Likelihood Estimation (MLE).

In short, that involves either:

1. using an analytical closed-form solution to fit the parameters exactly
2. solving an optimisation problem to maximise a log-likelihood function, or (more usually) minimise a negative log-likelihood function

count.data <- data.frame(events=0:13, count=c(26069, 30175, 18997, 8136, 2934, 820, 250, 54, 16, 4, 7, 1, 0, 1))


And we'll convert this count data to actual raw observations, for use later:

raw.data <- rep(count.data$events, count.data$count)


Before going on to the Negative Binomial, let's first fit a Poisson distribution to your data using MLE. This is actually very simple, because the ML estimate of the Poisson parameter $\lambda$ is simply the sample mean:

poisson.lambda.mle <- mean(raw.data)


We can use this value for $\lambda$ to generate Poisson counts, to compare against yours later:

fitted.poisson.counts <- dpois(count.data$events, poisson.lambda.mle) * sum(count.data$count)


The Negative Binomial distribution is much more complicated to fit. One reason is that there's no simple analytical expression for its parameters, so it needs to be formulated and solved as an optimisation problem. More importantly, if you parameterise it by an integer $r$, like you are doing, this integer parameter makes the optimisation much harder in R. In fact, I don't know of any native R integer programming functions.

So let's have a quick review of the Negative Binomial: there are actually quite a few ways to parameterise this distribution. The one you've chosen is a "coin-flipping" interpretation based on the number of success/failures, similar to the way the Binomial distribution is taught. There are a few such parameterisations but we actually won't use any of them, for the reasons already discussed.

Instead, the Negative Binomial distribution can be interpreted as a Gamma-mixture of Poisson distributions. This can also be parameterised in a few ways. We will use this one (as described by the R documentation):

An alternative parametrization (often used in ecology) is by the mean mu (see above), and size, the dispersion parameter, where prob = size/(size+mu). The variance is mu + mu^2/size in this parametrization.

The reason we will use this interpreation & parameterisation is because it allows us to use the MASS package's readymade implementation of MLE for the Negative Binomial. Let's get to it:

library(MASS)
negbin.fitted.object <- MASS::fitdistr(raw.data, "Negative Binomial")


Once again, we can use these ML estimates for the size and mu parameters to generate counts to compare against yours:

negbin.size.mle <- negbin.fitted.object$estimate['size'] negbin.mu.mle <- negbin.fitted.object$estimate['mu']
fitted.negbin.counts <- dnbinom(count.data$events, size=negbin.size.mle, mu=negbin.mu.mle) * sum(count.data$count)


And that's it. Not only have we fitted Poisson and Negative Binomial distributions to your data, we have prepared counts to compare against yours. Let's print them out together:

count.comparison <- transform(count.data,
fitted.poisson.counts=fitted.poisson.counts, fitted.negbin.counts=fitted.negbin.counts)
print(round(count.comparison, 0), row.names=FALSE)


which outputs:

 events count fitted.poisson.counts fitted.negbin.counts
0 26069                 24712                26149
1 30175                 31234                30178
2 18997                 19739                18727
3  8136                  8316                 8291
4  2934                  2628                 2933
5   820                   664                  881
6   250                   140                  233
7    54                    25                   56
8    16                     4                   12
9     4                     1                    3
10     7                     0                    0
11     1                     0                    0
12     0                     0                    0
13     1                     0                    0


As you rightly suspected, the counts generated by the ML Negative Binomial distribution are closer to your original count data than those of the ML Poisson distribution we have fitted.

This is even more apparent in a chart:

plot(count ~ events, count.data, 'b', col='black', main='Comparing actual & fitted counts')
lines(count.data$events, fitted.poisson.counts, 'b', col='red') lines(count.data$events, fitted.negbin.counts, 'b', col='blue')
grid()
legend("topright", legend=c('Actual counts', 'Fitted Poisson counts', 'Fitted Negative Binomial counts'), lty=1, col=c('black', 'red', 'blue'))


• Why has this been downvoted? Thanks very much for your answer, I am massively busy at the moment, but I will spend some time on analyzing your method over the next day or two. Thanks for your time. – MoonKnight Feb 13 '18 at 15:43
• No problem. I'd love to know why someone downvoted too. My answer does not explicitly formulate a log-likelihood and solve it, but the question didn't ask for that: it simply wants an R script to do the estimation, which is what my answer provides, together with a full explanation. – A. G. Feb 13 '18 at 15:47
• Perhaps a mistaken click from a cell phone review? One can only speculate... – jbowman Feb 13 '18 at 16:26
• @A.G. Sorry I have just gotten around to going through this - this is an AWESOME answer! Thanks very much for your help. But I wondered if I could ask you a few questions - I am happy to open up another question with bounty if required. The first question is (I have checked the MASS package documentation) what are the size and mu parameters in relation to the negative binomial distribution which usually uses r and p? – MoonKnight Apr 11 '18 at 20:42