2
$\begingroup$

I'm trying to fit a dataset to a Poisson distribution, but have probably messed up the parameters somewhere along the way. My data consists of 112 10 minute intervals where radiation hits a detector and is counted. Each 10 minute interval got ~1000 counts.

The issue is that after using scipy.optimize's curve_fit, I get essentially null values for all x (see picture).

import pandas as pd
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)
from scipy.special import factorial

############# This block is for loading data - the output, coinc, is length 112 array; [1125,1117,1056,...,1076]
df = pd.read_csv('CSMHUNT_2021-11-22_Overnight.csv',sep=';')
coinc=np.array(df[' COINC ' ][1:])
####################################################################################

def Normal(x,lamb):
    return norm.pdf(x,lamb,np.sqrt(lamb))

def Poisson(k,lamb):
    return np.exp(-lamb)*lamb**k *(1/factorial(k))

#histogram
plt.figure(dpi=100)
hist,bin_edges,params=plt.hist(coinc,bins=10,density=True,label='data')  #counts and bins

#Centres of bins
centers=(bin_edges[:-1]+bin_edges[1:])/2

#Guess for the lambda parameter. This is just the mean.
guess=np.mean(coinc)

par,cov = curve_fit(Poisson,centers,hist,p0=guess)
plt.plot(centers,Poisson(centers,*par),'r--',label='Fit')
plt.legend()

Output plot, with fitted function zero at all values.

I have a suspicion that I've gotten things turned around in my head, as the fit is obviously wrong somehow, but I can't spot the error. I tried replacing the starting guess lambda=np.mean(coinc) with np.mean(hist), which produces identically zero results.

p.s. I am aware that for for large lambda $(\lambda>1000)$ the normal distribution with mean $\lambda$ and variance $\lambda$ is a good approximation - hence why it's also defined in my code. It amounts to the same thing, however, as fitting to a Normal distribution still provides essentially null values (e-200 or so) for all the probabilities.

$\endgroup$
3
  • 3
    $\begingroup$ Why are you using curve fitting at all? What problem does it solve? The MLE for a Poisson is the mean, which you’ve already computed. $\endgroup$
    – Sycorax
    Nov 28 '21 at 19:01
  • $\begingroup$ I don't know python, but is it possible that curve_fit doesn't know that it should only evaluate the function at integers? (else factorial will give strange results) $\endgroup$
    – JDL
    Nov 29 '21 at 9:12
  • $\begingroup$ This histogram departs visibly from a Poisson shape (which, which this many counts, will be almost indistinguishable from a Normal distribution with a standard deviation around $33$ or so). Yes, you do have a programming problem, but fixing that isn't on topic here. $\endgroup$
    – whuber
    Nov 29 '21 at 14:29
6
$\begingroup$

I have two, non-exclusive hypotheses for the behavior.

  1. Floating point arithmetic is not sufficiently precise to represent large exponents and large factorials, causing catastrophic loss of precision.
  2. curve_fit isn't estimating the quantity that you want. There's simply no need to use the curve_fit function for this problem, because Poisson MLEs are easily computed.

This is fine, since we can just use the scipy functions for the Poisson distribution. The MLE of the Poisson parameter is the sample mean. I don't have access to your data, so I can't compute the mean, but you've already computed the mean in any event.

from scipy.stats import poisson
import numpy as np
from scipy.special import gammaln

k_vals = np.array([975, 1000, 1025, 1050, 1075, 1100,1125, 1150])
guess = 1030.4
probs = poisson.pmf(k_vals, mu=guess)

This gives reasonable results.

2.80323690e-03 8.01994963e-03 1.22844824e-02 1.02289554e-02
 4.69781334e-03 1.20658551e-03 1.75611575e-04 1.46676596e-05

Alternatively, we can write a quick-and-dirty log-scale implementation of the Poisson pmf and then exponentiate.

def dirty_poisson_pmf(x, mu):
    out = -mu + x * np.log(mu) - gammaln(x + 1)
    return np.exp(out)

dirty_probs = dirty_poisson_pmf(k_vals, mu=guess)
diff = probs - dirty_probs

And the differences are all on the order of machine epsilon for double-precision floating point.

[ 1.27502175e-15 -3.64638875e-15  5.58580959e-15  4.65252836e-15
 -2.13631196e-15  5.48606299e-16 -7.98514900e-17  6.66953743e-18]

When I try to directly apply the formula for the Poisson PMF, I get RuntimeWarnings from numpy that there has been overflow. The results are NAN. So we know that hypothesis (1) is correct.

Because hypothesis (2) is not necessary to estimate a Poisson PMF, and we've already, we don't need to worry about debugging curve_fit.

$\endgroup$
2
  • 2
    $\begingroup$ You import gammaln -- that's the trick! Work on the log scale whenever you have extremely large or small positive numbers. $\endgroup$
    – Sycorax
    Nov 28 '21 at 19:48
  • 1
    $\begingroup$ I deleted my previous comment, because I believe I've found the (actual) issue; the "centers" list in the above is NOT a list of integers - so the factorial function defaults to zero. If I instead use something like np.arange(int(min(centers)),int(max(centers))) it is computed just fine. $\endgroup$ Nov 28 '21 at 20:06

Not the answer you're looking for? Browse other questions tagged or ask your own question.