Fitting un-normalized gaussian in histogram python

Question

I have a dark image (raw format), and plotted the image and distribution of the image. As you can see, there is a peak at 16, please ignore that. I want to fit at gaussian curve through this histogram. I've used this method to fit: Un-normalized Gaussian curve on histogram. However; my Gaussian fit never comes close to what it supposed to be. Am I doing something wrong with turning the image into the correct format for the plot or is something else going wrong?

This is the current code I use to generate this data:

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

def fitGaussian(x,a,mean,sigma):
    return (a*np.exp(-((x-mean)**2/(2*sigma))))

fname = 'filepath.raw'
im = np.fromfile(fname,np.int16)
im.resize([3056,4064])

plt.figure()
plt.set_cmap(viridis)
plt.imshow(im, interpolation='none', vmin=16, vmax=np.percentile(im.ravel(),99))
plt.colorbar()
print 'Saving: ' + fname[:-4] + '.pdf'
plt.savefig(fname[:-4]+'.pdf')

plt.figure()
data = plt.hist(im.ravel(), bins=4096, range=(0,4095))

x = [0.5 * (data[1][i] + data[1][i+1]) for i in xrange(len(data[1])-1)]
y = data[0]

popt, pcov = curve_fit(fitGaussian, x, y, [500000,80,10])
x_fit = py.linspace(x[0], x[-1], 1000)
y_fit = fitGaussian(x_fit, *popt)
plt.plot(x_fit, y_fit, lw=4, color="r")        

plt.xlim(0,300)
plt.ylim(0,1e6)
plt.show()

EDIT: (response to Reblochon Masque)

If I remove the bin at 16 I still get the same fit:

I am not entirely sure I understand what you are doing, but it seems that for a fit you may want to eliminate the spike around 20 on the histogram — Reblochon Masque
@ReblochonMasque Even when I eliminate the bin at 16. I get a similar result, it only shifted to the right. See my edit. — SjonTeflon
I guess you get this curve because you fit the Gaußian to all the bins, most of which are zero. You could fit the curve to the non-zero bins only, but this is messy. — MB-F
It is already better... Good suggestion kazemakase; an easy way to check it this could be the case is to reduce your number of bins so all have values (try with 5-8 maybe?) — Reblochon Masque

MB-F MB-F · Accepted Answer · 2016-04-06T08:09:22

The fitted Gaussian appears too low because it is fit to all the bins, most of which are zero. A solution is to fit the Gaussian only to the non-zero bins.

I use np.histogram instead of plt.hist to get the bin vaules, but this is just a matter of taste. The important part is the definition of xh and yh

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

# Generate example data
x = np.random.randn(100000) * 50 + 75
x = np.round(x / 10) * 10
x = x[x >= 20]

yhist, xhist = np.histogram(x, bins=np.arange(4096))

xh = np.where(yhist > 0)[0]
yh = yhist[xh]

def gaussian(x, a, mean, sigma):
    return a * np.exp(-((x - mean)**2 / (2 * sigma**2)))

popt, pcov = curve_fit(gaussian, xh, yh, [10000, 100, 10])

plt.plot(yhist)
i = np.linspace(0, 300, 301)
plt.plot(i, gaussian(i, *popt))
plt.xlim(0, 300)

P.S. Sigma usually denotes the standard deviation rather than the variance. That is why I squared it in the gaussian function.

Fitting un-normalized gaussian in histogram python

1 Answers