У меня проблема с поиском метода наименьших квадратов для набора заданных данных. Я знаю, что данные следуют функции, ведьма представляет собой свертку гаусса и прямоугольника (рентгеновский снимок через широкую щель). Что я сделал до сих пор, так это взглянул на интеграл свертки и обнаружил, что он сводится к следующему: 
параметр интегрирования a - это ширина щели (неизвестная и желаемая) с g(xt) функцией Гаусса, определенной как 
Таким образом, в основном функция, которая подходит, представляет собой интеграционную функцию Гаусса с границами интегрирования, заданными параметром ширины a. Затем интегрирование также выполняется со сдвигом x-t.
Вот меньшая часть данных и ручная посадка. из pylab import * из scipy.optimize import curve_fit из scipy.integrate import 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
Чтобы убедиться, что 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 функции, но здесь он, похоже, не работает. С другой стороны, возможно, здесь я переоцениваю силу аппроксимации наименьших квадратов, и в данном случае она может оказаться непригодной, но мне не нравится использовать здесь метод максимального правдоподобия. В общем, я никогда не пытался подогнать интегральную функцию к данным, так что это для меня ново. Скорее всего проблема здесь. Мои знания о квадроциклах ограничены, и может быть лучший способ. Выполнение интеграла аналитически, насколько мне известно, невозможно, но явно было бы идеальным решением;).
Итак, есть идеи, откуда эта ошибка?
