I have got a frequency table of how many events occur within a 5-minute time window.

plotdata <- structure(list(events = 0:6, 
                           N = c(511468L, 75194L, 7813L, 1102L, 174L, 86L, 23L)), 
                      row.names = c(NA, -7L), class = c("data.table", "data.frame"))

#   events frequency
#1:      0    511468
#2:      1     75194
#3:      2      7813
#4:      3      1102
#5:      4       174
#6:      5        86
#7:      6        23

Explained: 7813 times there was a 5-minute timewindow in which two events occured.

I am trying to fit this data onto a poisson-curve, but I'm geting lost here. My last statistics course was waaay back in college and I'm getting lost in the terminology.

What I've tried so far:

find lambda:

lambda <- sum( plotdata$events * plotdata$N ) / sum( plotdata$N )
#lambda = 0.1600879401

Get estimated values of poisson-distribution Resulting in (too?) low estimated values for events > 2

plotdata[, poisson.P := exp( -1 * lambda ) * lambda^events / factorial( events )]
plotdata[, poisson.N := poisson.P * sum( N ) ][]

ggplot( plotdata, aes( x = events ) ) + 
  geom_line( aes( y = N ) ) + 
  geom_line( aes( y = poisson.N), colour = "red" ) + 

enter image description here

black = counted values, red = result from poisson

Am I doing something wrong here? Or is my data not suited for a description by poisson-distribution, or..., or... ?

Underestimation on a larger number of events is a no-go in my usecase. So I would really like the estimated output to perform better on events > 3

  • 1
    $\begingroup$ Your data are indeed not like a Poisson distribution; they're considerably more like a geometric. For example, a Poissonness plot shows a distinct "kink". You can see in your plot that the black line is somewhat nearer a straight line than it is like the red curve. $\endgroup$ – Glen_b -Reinstate Monica Sep 24 at 10:23
  • $\begingroup$ Could you elaborate on what you need to estimate? In particular, what would be the limitation of using the frequencies you have observed in this large dataset? $\endgroup$ – whuber Sep 24 at 11:38
  • $\begingroup$ @whuber I would like to be able to estimate the probability of events >6 (since they did not occur dring measuring period). $\endgroup$ – Wimpel Sep 24 at 13:50
  • $\begingroup$ That's going to be problematic, because you are basically having to guess what that probability might be: you have no data (apart from the fact that you haven't recorded any such values in a dataset of almost 600K observations). Anything you do will necessarily be an extrapolation beyond the data you have and therefore will depend strongly on your assumptions about how the probabilities change. $\endgroup$ – whuber Sep 24 at 13:56
  • $\begingroup$ @whuber, ahah.. i see... I was hopig to fit a (poisson) distribution on the observed data, and then calculate odds for unobserved data. Also, if I can find the correct distribution, I can (probably) use it into a MonteCarlo analysis later on in my project. $\endgroup$ – Wimpel Sep 24 at 14:10

This data do not seem to be following a Poisson distribution, and to me it is not clear that zero inflation is the issue. The negative binomial does give a better fit. One way to see this is (following your R code in the post):

 (lambda <- sum( plotdata$events * plotdata$N ) / sum( plotdata$N ) )
[1] 0.1600879
> (var    <- Hmisc::wtd.var(plotdata$events, plotdata$N) )  
[1] 0.1793297

A visualization for Poissonness:

vcd::distplot(cbind(plotdata$N, plotdata$events), "poisson") 

Poissonness plot

And ditto for the negative binomial:

vcd::distplot(cbind(plotdata$N, plotdata$events), "nbinomial")

negbinness plot

but the fit is better, not perfect.

  • 1
    $\begingroup$ Thanks! this gives me a nice new starting point to work with! Time to hit the books/tutorials again on the binominal distribution :) $\endgroup$ – Wimpel Sep 24 at 18:02

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.