Solved – Fitting pmf of a scaled Poisson distribution and Python histogram plotting

Question

I have a nuclei meanlife of $550\mu s$, for which I've taken the frequency(rate) to be $1/meanlife = 1818$. I then sampled randomly from a poisson distribution with that frequency, taking the reciprocal of the sample and plotted it on a histogram. Given that the frequency will be distributed poissonly according to $$P(x) = \frac{e^{-\lambda}\lambda^x}{x!}$$
Since I'm plotting the histogram of $t = 1/x$ where I'm sampling x randomly from a Poisson distribution, I thought I'd fit a line of
$$P(t) = \frac{e^{-\lambda}\lambda^{\frac{1}{t}}}{\frac{1}{t}!}$$
Not sure if the Jacobian is needed but if it was then it'd just be a matter of multiplying $P(t)$ by $1/t^2$.
The whole code in python looks something like this

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.special import factorial
meanlife = 550e-6
decay_lifetimes = 1/np.random.poisson((1/meanlife), size=100000) # I'm taking 100000 independent samples basically

def poisson(t, rate): #x axis is k
    return (rate**(1/t)/factorial(1/t))*np.exp(-rate)

hist, bins = np.histogram(decay_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.0001, [5e-4]))
plt.plot(center, poisson(center, *popt), 'r--', label='Poisson fit')
plt.legend(loc = 'best')
plt.tight_layout()

which gives me a histogram like this.

I get the correct histogram which is what I expected. However as you can see the fitted red line on it keeps giving me 0 everywhere, because of my frequency(rate) of 1818, which suggests that my transformation of pmf is wrong. Not sure what I could do here to remedy the problem.

Answer 1

In your example the rate is large (>1000) and in this case the normal distribution with mean $\lambda$, variance $\lambda$ is a very good approximation to the poisson with rate $\lambda$. So you could consider fitting a normal to your data instead. It is also important to choose an appropriate initial value for the parameter. (Otherwise, the default initial value is 1, which is not a very good guess for your data.)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import curve_fit
from scipy.special import gammaln # x! = Gamma(x+1)

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

def transformation_and_jacobian(x):
    return 1./x, 1./x**2.

def tfm_normal_pdf(x, lam):
    y, J = transformation_and_jacobian(x)
    return norm.pdf(y, lam, np.sqrt(lam)) * J

def tfm_poisson_pdf(x, mu):
    y, J = transformation_and_jacobian(x)
    # For numerical stability, compute exp(log(f(x)))
    return np.exp(y * np.log(mu) - mu - gammaln(y + 1.)) * J

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

# Important: Choose a reasonable starting point
p0 = 1 / np.mean(decay_lifetimes)

norm_opt, _ = curve_fit(tfm_normal_pdf, center, hist, p0=p0)
pois_opt, _ = curve_fit(tfm_poisson_pdf, center, hist, p0=p0)

plt.plot(center, tfm_normal_pdf(center, *norm_opt), 'g--', label='(Transformed) Normal fit')
plt.plot(center, tfm_poisson_pdf(center, *pois_opt), 'r--', label='(Transformed) Poisson fit')
plt.legend(loc = 'best')
plt.tight_layout()