在直方图中拟合非归一化高斯 python
Fitting un-normalized gaussian in histogram python
我有一个深色图像(原始格式),并绘制了图像和图像的分布。如您所见,在 16 处有一个峰值,请忽略它。我想通过这个直方图拟合高斯曲线。我用这种方法来适应:
Un-normalized Gaussian curve on histogram。然而;我的高斯拟合永远不会接近它应该的样子。将图像转换为正确的绘图格式是我做错了什么,还是其他地方出了问题?
这是我用来生成此数据的当前代码:
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()
编辑:(对 Reblochon 面膜的回应)
如果我在 16 点移除垃圾桶,我仍然得到同样的合身:
拟合的高斯分布看起来太低,因为它适合所有 bin,其中大部分为零。一种解决方案是仅将高斯拟合到非零区间。
我使用 np.histogram
而不是 plt.hist
来获取 bin vaules,但这只是个人喜好问题。重要的部分是 xh
和 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 通常表示标准偏差而不是方差。这就是为什么我在 gaussian
函数中对它求平方。
我有一个深色图像(原始格式),并绘制了图像和图像的分布。如您所见,在 16 处有一个峰值,请忽略它。我想通过这个直方图拟合高斯曲线。我用这种方法来适应:
Un-normalized Gaussian curve on histogram。然而;我的高斯拟合永远不会接近它应该的样子。将图像转换为正确的绘图格式是我做错了什么,还是其他地方出了问题?
这是我用来生成此数据的当前代码:
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()
编辑:(对 Reblochon 面膜的回应)
如果我在 16 点移除垃圾桶,我仍然得到同样的合身:
拟合的高斯分布看起来太低,因为它适合所有 bin,其中大部分为零。一种解决方案是仅将高斯拟合到非零区间。
我使用 np.histogram
而不是 plt.hist
来获取 bin vaules,但这只是个人喜好问题。重要的部分是 xh
和 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 通常表示标准偏差而不是方差。这就是为什么我在 gaussian
函数中对它求平方。