Fitting to Poisson histogram

Question

I am trying to fit a curve over the histogram of a Poisson distribution that looks like this

I have modified the fit function so that it resembles a Poisson distribution, with the parameter t as a variable. But the curve_fit function can not be plotted and I am not sure why.

def histo(bsize):
    N = bsize
    #binwidth
    bw = (dt.max()-dt.min())/(N-1.)
    bin1 = dt.min()+ bw*np.arange(N)
    #define the array to hold the occurrence count
    bincount= np.array([])
    for bin in bin1:
        count = np.where((dt>=bin)&(dt<bin+bw))[0].size
        bincount = np.append(bincount,count)
    #bin center
    binc = bin1+0.5*bw
    plt.figure()
    plt.plot(binc,bincount,drawstyle= 'steps-mid')
    plt.xlabel("Interval[ticks]")
    plt.ylabel("Frequency")
histo(30)
plt.xlim(0,.5e8)
plt.ylim(0,25000)
import numpy as np
from scipy.optimize import curve_fit
delta_t = 1.42e7
def func(x, t):
    return t * np.exp(- delta_t/t) 
popt, pcov = curve_fit(func, np.arange(0,.5e8),histo(30))
plt.plot(popt)

Answer 1

The problem with your code is that you do not know what the return values of curve_fit are. It is the parameters for the fit-function and their covariance matrix - not something you can plot directly.

Binned Least-Squares Fit

In general you can get everything much, much more easily:

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit from scipy.special import factorial from scipy.stats import poisson  # get poisson deviated random numbers data = np.random.poisson(2, 1000)  # the bins should be of integer width, because poisson is an integer distribution bins = np.arange(11) - 0.5 entries, bin_edges, patches = plt.hist(data, bins=bins, density=True, label='Data')  # calculate bin centres bin_middles = 0.5 * (bin_edges[1:] + bin_edges[:-1])   def fit_function(k, lamb):     '''poisson function, parameter lamb is the fit parameter'''     return poisson.pmf(k, lamb)   # fit with curve_fit parameters, cov_matrix = curve_fit(fit_function, bin_middles, entries)  # plot poisson-deviation with fitted parameter x_plot = np.arange(0, 15)  plt.plot(     x_plot,     fit_function(x_plot, *parameters),     marker='o', linestyle='',     label='Fit result', ) plt.legend() plt.show() 

This is the result:

Unbinned Maximum-Likelihood fit

An even better possibility would be to not use a histogram at all and instead to carry out a maximum-likelihood fit.

But by closer examination even this is unnecessary, because the maximum-likelihood estimator for the parameter of the poissonian distribution is the arithmetic mean.

However, if you have other, more complicated PDFs, you can use this as example:

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize from scipy.special import factorial from scipy import stats   def poisson(k, lamb):     """poisson pdf, parameter lamb is the fit parameter"""     return (lamb**k/factorial(k)) * np.exp(-lamb)   def negative_log_likelihood(params, data):     """     The negative log-Likelihood-Function     """      lnl = - np.sum(np.log(poisson(data, params[0])))     return lnl  def negative_log_likelihood(params, data):     ''' better alternative using scipy '''     return -stats.poisson.logpmf(data, params[0]).sum()   # get poisson deviated random numbers data = np.random.poisson(2, 1000)  # minimize the negative log-Likelihood  result = minimize(negative_log_likelihood,  # function to minimize                   x0=np.ones(1),            # start value                   args=(data,),             # additional arguments for function                   method='Powell',          # minimization method, see docs                   ) # result is a scipy optimize result object, the fit parameters  # are stored in result.x print(result)  # plot poisson-distribution with fitted parameter x_plot = np.arange(0, 15)  plt.plot(     x_plot,     stats.poisson.pmf(x_plot, *parameters),     marker='o', linestyle='',     label='Fit result', ) plt.legend() plt.show()