CODEWITHSUNDEEP
python-3.xstatisticsdistributioncurve-fittingpoisson
user3613025
user3613025

Reputation: 383

Fitting Poisson distribution on a histogram

Despite the overwhelming amount of posts on fitting Poisson distribution onto a histogram, having followed all of them, none of them seems to work for me.

I'm looking to fit a poisson distribution on this histogram which I've plotted as such:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.misc import factorial

def poisson(t, rate, scale): #scale is added here so the y-axis 
# of the fit fits the height of histogram
    return (scale*(rate**t/factorial(t))*np.exp(-rate))

lifetimes = 1/np.random.poisson((1/550e-6), size=100000)

hist, bins = np.histogram(lifetimes, bins=50)
width = 0.8*(bins[1]-bins[0])
center = (bins[:-1]+bins[1:])/2
plt.bar(center, hist, align='center', width=width, label = 'Normalised data')

popt, pcov = curve_fit(poisson, center, hist, bounds=(0.001, [2000, 7000]))
plt.plot(center, poisson(center, *popt), 'r--', label='Poisson fit')
# import pdb; pdb.set_trace()
plt.legend(loc = 'best')
plt.tight_layout()

The histogram I get looks like this:

enter image description here

I gave the guess of scale as 7000 to scale the distribution to the same height as the y-axis of the histogram I plotted and a guess of 2000 as the rate parameter since it's 2000 > 1/550e-6. As you can see the fitted red dotted line is 0 at every point. Weirdly pdb.set_trace() tells me that the poisson(center, *popt) gives me a list of 0 values.

    126     plt.plot(center, poisson(center, *popt), 'r--', label='Poisson fit')
    127     import pdb; pdb.set_trace()
--> 128     plt.legend(loc = 'best')
    129     plt.tight_layout()
    130 


ipdb> 
ipdb> poisson(center, *popt)
array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])

Which doesn't make sense. What I want is to fit a poisson distribution on the histogram such that it finds the best coefficient of the poisson distribution equation. I suspected that it might be have to do with because I am plotting histogram of lifetimes instead, which is technically randomly sampled data from the inverse of the poisson distribution. So I tried to compute the jacobian of the distribution so I can make a change of variables but it still won't work. I feel like I'm missing something here that's not coding but rather mathematics related.

Upvotes: 1

Views: 1117

Answers (1)

Nathan Hamm
Nathan Hamm

Reputation: 300

You're calculation is rounding to zero. With a rate of 2000 and scale of 7000 your poisson formula is reduced to:

7000 * 2000^t/(e^(2000) * t!)

Using Stirling's approximation t! ~ (2*pi*t)^(1/2) * (t/e)^t you get:

[7000 * 2000^t] / [Sqrt(2*pi*t) * e^(2000-t) * (t^t)] ~ poisson(t)

I used python to get the first couple values of poisson(t):

poisson(1) -> 0
poisson(2) -> 0
poisson(3) -> 0

Using wolfram alpha you find that the derivative of the denominator greater than the derivative of the numerator for all real numbers greater than zero. Therefore, poisson(t) is approaching zero as t gets larger.

This means that no matter what t is, if your rate is 2000, the poisson function will return 0.

Sorry for the formatting. They wont let me post TeX yet.

Upvotes: 0

Related Questions