CODEWITHSUNDEEP
pythonmatplotlibcurve-fittingdata-fitting
Yababaa
Yababaa

Reputation: 577

Fitting un-normalized gaussian in histogram python

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? enter image description here Gaussian distribution of the image

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: enter image description here

Upvotes: 3

Views: 7402

Answers (1)

MB-F
MB-F

Reputation: 23637

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)

enter image description here

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

Upvotes: 5

Related Questions