首页 文章

在直方图python中拟合非标准化高斯

提问于
浏览
2

我有一个暗图像(原始格式),并绘制图像的图像和分布 . 正如您所看到的那样,在16处有一个高峰,请忽略它 . 我想通过这个直方图拟合高斯曲线 . 我已经使用这种方法来适应:Un-normalized Gaussian curve on histogram . 然而;我的高斯拟合永远不会接近它应该是什么 . 将图像转换为正确的图形格式或者出现其他问题时,我做错了吗?
enter image description here

Gaussian distribution of the image

这是我用来生成此数据的当前代码:

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: (对Reblochon Masque的回应)

如果我在16处移除垃圾箱,我仍然可以获得相同的适合度:
enter image description here

1 回答

  • 3

    拟合高斯看起来太低,因为它适合所有的箱子,其中大多数是零 . 一种解决方案是仅将高斯拟合到非零二进制位 .

    我使用 np.histogram 而不是 plt.hist 来获得bin vaules,但这只是一个品味问题 . 重要的部分是 xhyh 的定义

    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

    附: Sigma通常表示标准偏差而不是方差 . 这就是我在 gaussian 函数中将其平方的原因 .

相关问题