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)][1].

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

Let's first load up your count data in a dataframe:

    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][2]:

    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][3] 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][4] 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][5]. This can also be parameterised in a few ways. We will use this one ([as described by the R documentation][6]):

> 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][7]. 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'))

[![enter image description here][8]][8]


  [1]: https://en.wikipedia.org/wiki/Maximum_likelihood_estimation
  [2]: https://www.statlect.com/fundamentals-of-statistics/Poisson-distribution-maximum-likelihood
  [3]: https://en.wikipedia.org/wiki/Integer_programming
  [4]: https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations
  [5]: https://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma%E2%80%93Poisson_mixture
  [6]: https://stats.stackexchange.com/questions/325625/establishing-the-probability-distribution-that-governs-a-random-process
  [7]: https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/fitdistr.html
  [8]: https://i.sstatic.net/8HHaM.png