如何根据数据找到合适的高斯曲线？

Question

我有一段代码可以从一维高斯分布中生成数据：

import numpy as np
from scipy.stats import norm
from scipy.optimize import curve_fit
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import gauss

# Beginning in one dimension:
mean = 0; Var = 1; N = 1000
scatter = np.random.normal(mean,np.sqrt(Var),N)
scatter = np.sort(scatter)
mu,sigma = norm.fit(scatter)

我通过 norm.fit() 来获取均值（mu）和标准差（sigma）。

现在我想用

xdata = np.linspace(-5,5,N)
pop, pcov = curve_fit(gauss.gauss_1d,xdata,scatter)

来获取我的参数。

问题是，我不知道怎么把我那些散乱的点（从一维高斯分布中生成的）映射到 x 轴上，以便使用 curve_fit。

另外，假设我还是像之前那样直接使用 mu 和 sigma。

我用以下代码来绘图：

n, bins, patches = plt.hist(scatter,50,facecolor='green')
y = 2*max(n)*mlab.normpdf(bins,mu,sigma)
l = plt.plot(bins,y,'r--')

plt.xlabel('x-coord')
plt.ylabel('Occurrences')
plt.grid(True)
plt.show()

但我得猜测振幅（amplitude），用的是 2*max(n)。这样做虽然能用，但不太可靠。我该怎么找到振幅，而不需要猜测呢？

Answer 1

为了避免猜测幅度，可以在调用 hist() 时设置 normed=True，这样得到的幅度就对应于 normpdf() 的值。

如果要进行曲线拟合，我建议使用累积分布而不是密度：每个样本的高度是 1/N，这些高度会逐渐加起来到1。这样做的好处是你不需要把样本分组到不同的区间里。

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

# Beginning in one dimension:
mean = 0; Var = 1; N = 100
scatter = np.random.normal(mean,np.sqrt(Var),N)
scatter = np.sort(scatter)
mu1,sigma1 = norm.fit(scatter) # classical fit

scat_sum = np.cumsum(np.ones(scatter.shape))/N # cumulative samples
[mu2,sigma2],Cx = curve_fit(norm.cdf, scatter, scat_sum, p0=[0,1]) # curve fit
print(u"norm.fit():  µ1= {:+.4f}, σ1={:.4f}".format(mu1, sigma1))
print(u"curve_fit(): µ2= {:+.4f}, σ2={:.4f}".format(mu2, sigma2))

fg = plt.figure(1); fg.clf()
ax = fg.add_subplot(1, 1, 1)
t = np.linspace(-4,4, 1000)
ax.plot(t, norm.cdf(t, mu1, sigma1), alpha=.5, label="norm.fit()")
ax.plot(t, norm.cdf(t, mu2, sigma2), alpha=.5, label="curve_fit()")
ax.step(scatter, scat_sum, 'x-', where='post', alpha=.5, label="Samples")
ax.legend(loc="best")
ax.grid(True)
ax.set_xlabel("$x$")
ax.set_ylabel("Cumulative Probability Density")
ax.set_title("Fit to Normal Distribution")

fg.canvas.draw()
plt.show()

打印结果

norm.fit():  µ1= +0.1534, σ1=1.0203
curve_fit(): µ2= +0.1135, σ2=1.0444

并绘制图形