我在寻找一组给定数据的最小二乘法时遇到了一个问题。 我知道数据遵循一个函数,它是高斯和矩形的卷积(x射线穿过宽狭缝)。到目前为止,我所做的是看一看卷积积分,发现它可以归结为: 积分参数a是狭缝的宽度(未知和期望的),g(x-t)高斯函数定义为 所以基本上拟合的函数是一个高斯积分函数,积分边界由宽度参数a给出。积分也在x-t偏移的情况下进行
这是一个较小的部分数据和手工拟合。 从pylab进口* 从scipy.optimize公司导入曲线拟合 从整合导入四边形
# 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()
^{pr2}$
为了证明func确实描述了数据,并且我的计算是正确的,我反复研究了数据和函数,并试图匹配它们。 我发现以下是可行的:
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中的向量函数作为积分函数(func<;>vfunc),但是不使用它也没有帮助。一般来说,可以将numpy数组传递给函数,但在这里它似乎不起作用。另一方面,我可能高估了leas平方拟合的能力,在这种情况下它可能不可用,但我不喜欢在这里使用最大似然法。总的来说,我从来没有尝试过让一个积分函数适合数据,所以这对我来说是新的。问题可能就在这里。我对quad的了解有限,可能还有更好的方法。据我所知,用分析的方法进行积分是不可能的,但显然是理想的解决方案;)。在
你知道这个错误是从哪里来的吗?在
你有两个问题。一个是
quad
返回一个包含值和错误估计值的元组,另一个是如何向量化。您不想对vectors参数进行矢量化。np.vectorize
有一个for循环,因此您自己这样做不会带来性能提升:但没有从cd3中拿走。另外,我选择
quad
的第一个输出。在虽然这解决了你的问题,为了适应卷积,你可以考虑去傅立叶空间。卷积的傅里叶变换是函数变换的乘积,这将大大简化你的生活。此外,一旦进入傅立叶空间,你可以考虑应用低通滤波器来降低噪声。210个数据点足够高,可以获得良好的结果。在
另外,如果您需要更强大的算法,您应该考虑iminuit,使用ROOT久经考验的细节。在
相关问题 更多 >
编程相关推荐