将高斯积分函数拟合到数据上
我在寻找一组给定数据的最小二乘拟合时遇到了问题。我知道这些数据遵循一个函数，这个函数是高斯函数和矩形函数的卷积（就像X射线通过一个宽缝隙）。到目前为止，我查看了卷积积分，发现它可以简化为这个样子：
这里的积分参数a是缝隙的宽度（这是未知的，也是我想要找的），而g(x-t)是一个定义好的高斯函数，表示为：
所以，基本上要拟合的函数是一个高斯函数的积分，积分的边界由宽度参数a给出。积分时还会有一个x-t的偏移。
这是数据的一小部分和我手动拟合的结果。
从pylab导入所有内容
从scipy.optimize导入curve_fit
从scipy.integrate导入quad
# 1/10 of the Data to show the form.
xData = array([-0.1 , -0.09, -0.08, -0.07, -0.06, -0.05, -0.04, -0.03, -0.02,
-0.01, 0. , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07,
0.08, 0.09, 0.1 ])
yData = array([ 18. , 22. , 22. , 34.000999, 54.002998,
152.022995, 398.15799 , 628.39502 , 884.781982, 848.719971,
854.72998 , 842.710022, 762.580994, 660.435974, 346.119995,
138.018997, 40.001999, 8. , 6. , 4. ,
6. ])
yerr = 0.1*yData # uncertainty of the data
plt.scatter(xData, yData)
plt.show()
# functions
def gaus(x, *p):
""" gaussian with p = A, mu, sigma """
A, mu, sigma = p
return A/(sqrt(2*pi)*sigma)*numpy.exp(-(x-mu)**2/(2.*sigma**2))
def func(x,*p):
""" Convolution of gaussian and rectangle is a gaussian integral.
Parameters: A, mu, sigma, a"""
A, mu, sigma, a = p
return quad(lambda t: gaus(x-t,A,mu,sigma),-a,a)
vfunc = vectorize(func) # Probably this is a Problem but if I dont use it, func can only be evaluated at 1 point not an array
为了验证这个函数确实描述了数据，并且我的计算是正确的，我尝试了不同的数据和函数，想要让它们匹配。我发现以下方法是可行的：
p0=[850,0,0.01, 0.04] # will be used as starting values for fitting
sample = linspace(-0.1,0.1,200) # just to make the plot smooth
y, dy = vfunc(sample,*p0)
plt.plot(sample, y, label="Handmade Fit")
plt.scatter(xData, yData, label="Data")
plt.legend()
plt.show()
当我尝试使用刚获得的初始值来拟合数据时，问题就出现了：
fp, Sfcov = curve_fit(vfunc, xData, yData, p0=p0, sigma=yerr)
yf = vfunc(xData, fp)
plt.plot(x, yf, label="Fit")
plt.show()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-83-6d362c4b9204> in <module>()
----> 1 fp, Sfcov = curve_fit(vfunc, xData, yData, p0=p0, sigma=yerr)
2 yf = vfunc(xData,fp)
3 plt.plot(x,yf, label="Fit")
/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, **kw)
531 # Remove full_output from kw, otherwise we're passing it in twice.
532 return_full = kw.pop('full_output', False)
--> 533 res = leastsq(func, p0, args=args, full_output=1, **kw)
534 (popt, pcov, infodict, errmsg, ier) = res
535
/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
369 m = shape[0]
370 if n > m:
--> 371 raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
372 if epsfcn is None:
373 epsfcn = finfo(dtype).eps
TypeError: Improper input: N=4 must not exceed M=2
我觉得这意味着我拥有的数据点比拟合参数要少。我们来看看：
print("Fit-Parameters: %i"%len(p0))
print("Datapoints: %i"%len(yData))
Fit-Parameters: 4
Datapoints: 21
实际上，我有210个数据点。
正如上面所说的，我不太明白为什么在积分函数中需要使用numpy的vectorise函数（func <> vfunc），而不使用它也没有帮助。一般来说，可以将numpy数组传递给一个函数，但在这里似乎不太奏效。另一方面，我可能高估了最小二乘拟合的能力，它可能在这种情况下不适用，但我不想使用最大似然法。总的来说，我从未尝试过将积分函数拟合到数据上，所以这对我来说是新的。问题可能出在这里。我对quad的了解有限，可能还有更好的方法。根据我的知识，解析地进行积分是不可能的，但显然这是理想的解决方案；）
所以，有没有人能告诉我这个错误是从哪里来的？
1 个回答
你遇到了两个问题。第一个是 quad 函数返回的是一个元组，里面包含了计算结果和一个误差估计；第二个问题是在你处理向量时的方式。你不应该在向量参数上进行向量化。因为 np.vectorize 其实是用一个循环来处理的，所以自己做并不会提高性能：
def func(x, p):
""" Convolution of gaussian and rectangle is a gaussian integral.
Parameters: A, mu, sigma, a"""
A, mu, sigma, a = p
return quad(lambda t: gaus(x-t,A,mu,sigma),-a,a)[0]
def vfunc(x, *p):
evaluations = numpy.array([func(i, p) for i in x])
return evaluations
注意，我把 func 中的 * 去掉了，但 gaus 中的没有去掉。同时，我选择了 quad 的第一个输出结果。
虽然这样解决了你的问题，但如果你想做卷积，建议你可以考虑转到傅里叶空间。卷积的傅里叶变换是函数变换的乘积，这样会让你的工作简单很多。此外，一旦进入傅里叶空间，你还可以考虑使用低通滤波器来减少噪声。210个数据点已经足够让你得到不错的结果。
另外，如果你需要更强大的算法，可以考虑使用 iminuit，它是基于 ROOT 的 Minuit，经过了长时间的验证。