Revisions to Establishing the Probability Distribution that Governs a Random Process

added 289 characters in body

Source Link

edited Jan 29, 2022 at 1:26

82.8k
32
201
663

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

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

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

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

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

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

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

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

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

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

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)

    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)

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)

    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)

 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

     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

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

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

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

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

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

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

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

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)

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)

 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

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

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

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

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

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

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

    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)

    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)

     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

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

Bounty Ended with 50 reputation awarded by CommunityBot

occurred Feb 16, 2018 at 21:35

makes code runnable: prefix `count.` was missing, now added.

Source Link

edit approved Feb 13, 2018 at 18:20

Jim

2.2k
2
17
20

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

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(data$$events, size=negbin.size.mle, mu=negbin.mu.mle) * sum(count.data$count)

plot(count ~ events, count.data, 'b', col='black', main='Comparing actual & fitted counts')
lines(count.data$events, fitted.poisson.counts, 'b', col='red')
lines(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'))

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

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(data$count)

plot(count ~ events, data, 'b', col='black', main='Comparing actual & fitted counts')
lines(data$events, fitted.poisson.counts, 'b', col='red')
lines(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'))

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

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)

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

Source Link

created Feb 13, 2018 at 15:17

A. G.

2.2k
10
17

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:

using an analytical closed-form solution to fit the parameters exactly
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:

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(data$events, poisson.lambda.mle) * sum(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(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, data, 'b', col='black', main='Comparing actual & fitted counts')
lines(data$events, fitted.poisson.counts, 'b', col='red')
lines(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'))

Stack Exchange Network

Return to Answer