python - 对数据点拟合更好的高斯分布？

Question

我正在尝试将高斯拟合到一组似乎遵循高斯分布的数据点。我已经检查了很多可能的方法来做到这一点，但我并不真正理解其中的大多数。然而，我找到了一个似乎有效的解决方案，但我得到的实际拟合结果看起来并不比我的数据点更像高斯。

这是我的代码:

import numpy as np
import matplotlib.pyplot as plt
from scipy import asarray as ar, exp, sqrt
from scipy.optimize import curve_fit


angles = [-8, -6, -4, -2, 0, 2, 4, 6, 8]
data = [99, 610, 1271, 1804, 1823, 1346, 635, 125, 24]
angles = ar(angles)
data = ar(data)

n = len(x)
mean = sum(data*angles)/n
sigma = sqrt(sum(data*(angles-mean)**2)/n)

def gaus(x,a,mu,sigma):
    return a*exp(-(x-mu)**2/(2*sigma**2))

popt,pcov = curve_fit(gaus,angles,data,p0=[0.18,mean,sigma])


fig = plt.figure()
plt.plot(angles, data, "ob", label = "Measured")
plt.plot(angles,gaus(angles,*popt),'r',label='Fit')
plt.xlim(-10, 10)
plt.ylim(0, 2000)
plt.xticks(angles)
plt.title("$^{137}$Cs Zero Point")
plt.xlabel("Angle [$^\circ$]")
plt.ylabel("662 keV-Photon Count")
plt.grid()
plt.legend()
plt.show()

这是它生成的输出:

如您所见，拟合并没有描述出良好且对称的“真实”高斯分布。 有什么方法可以获得“更好”的高斯分布或者这已经是最好的了吗？

非常感谢!

Answer 1

我认为这里有两件事:

seem to follow a gaussian distribution

→ 如果您认为数据呈正态分布，那么您处于统计和概率分布领域，并且可能想要创建 test看看他们是否同意特定的分布(正态分布或其他分布)。

<小时/>

并处理你的情节:

get a "better" gaussian plot

在代码中，您可以省略 curve_fit 中的第一个估计，并针对连续自变量绘制拟合曲线:

import numpy as np
import matplotlib.pyplot as plt
from scipy import asarray as ar, exp, sqrt
from scipy.optimize import curve_fit


angles = [-8, -6, -4, -2, 0, 2, 4, 6, 8]
data = [99, 610, 1271, 1804, 1823, 1346, 635, 125, 24]
angles = ar(angles)
data = ar(data)

n = len(data)  ## <---
mean = sum(data*angles)/n
sigma = sqrt(sum(data*(angles-mean)**2)/n)

def gaus(x,a,mu,sigma):
    return a*exp(-(x-mu)**2/(2*sigma**2))

popt,pcov = curve_fit(gaus,angles,data)#,p0=[0.18,mean,sigma])  ## <--- leave out the first estimation of the parameters
xx = np.linspace( -10, 10, 100 )  ## <--- calculate against a continuous variable

fig = plt.figure()
plt.plot(angles, data, "ob", label = "Measured")
plt.plot(xx,gaus(xx,*popt),'r',label='Fit')  ## <--- plot against the contious variable
plt.xlim(-10, 10)
plt.ylim(0, 2000)
plt.xticks(angles)
plt.title("$^{137}$Cs Zero Point")
plt.xlabel("Angle [$^\circ$]")
plt.ylabel("662 keV-Photon Count")
plt.grid()
plt.legend()
plt.savefig('normal.png')
plt.show()

<小时/>

在此示例中:

print( popt )

[  1.93154077e+03  -9.21486804e-01   3.26251063e+00]

请注意，参数的第一次估计与结果相差几个数量级:0.18 与 1931.15。