将高斯积分函数拟合到数据

Posted 2023-03-12

【中文标题】将高斯积分函数拟合到数据

【英文标题】：Fit gaussian integral function to data

【发布时间】：2014-08-05 11:36:36

【问题描述】：

我在为一组给定数据找到最小二乘拟合时遇到问题。 我知道数据遵循一个函数，它是高斯和矩形的卷积（通过宽缝的 X 射线）。到目前为止，我所做的是查看卷积积分并发现它归结为： 积分参数 a 是狭缝的宽度（未知和期望），g(x-t) 是一个高斯函数，定义为 所以基本上要拟合的函数是高斯的积分函数，积分边界由宽度参数 a 给出。然后也以 x-t 的位移进行积分。

这是数据的一小部分和手工合身。 从 pylab 导入 * 从 scipy.optimize 导入曲线拟合 从 scipy.integrate 导入四边形

# 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

看到 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 中的 vectorise 函数作为积分函数 (func vfunc) 但不使用它也无济于事。一般来说，可以将一个 numpy 数组传递给一个函数，但它似乎在这里不起作用。另一方面，我可能高估了 leas-square-fit 的力量，在这种情况下它可能不可用，但我不喜欢在这里使用最大似然。一般来说，我从未尝试将积分函数拟合到数据中，所以这对我来说是新的。问题可能就在这里。我对四边形的了解有限，可能有更好的方法。据我所知，以分析方式进行积分是不可能的，但显然是理想的解决方案；）。

那么有什么想法是这个错误的来源吗？

【参考方案1】：

你有两个问题。一个是quad 返回一个包含错误值和估计值的元组，另一个是向量化的方式。您不想对向量参数进行向量化。 np.vectorize 有一个 for 循环，所以自己做不会有性能提升：

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 个数据点足以获得良好的结果。

另外，如果您需要更强大的算法，您应该考虑使用 ROOT 久经考验的 Minuit。

【讨论】：

完美。非常感谢。关于 ROOT 和傅立叶空间，您可能是对的。下次我深入研究它。到目前为止，我对您的更正感到非常满意。