Проблема
Ваша основная проблема заключается в том, что у вас есть строго неопределенная проблема с подгонкой.Подумайте об этом так: у вас есть три неизвестных, но только одна точка данных.Это похоже на решение для x, y, z, когда у вас есть только одно уравнение.Поскольку высота вашего гауссиана может изменяться независимо от его ширины, существует бесконечно много распределений, все с различной шириной, которые будут удовлетворять ограничениям вашей подгонки.
Более конкретно, ваши параметры a и sigmaоба могут изменить максимальную высоту распределения, что является в значительной степени единственной вещью, которая имеет значение с точки зрения достижения хорошего соответствия (по крайней мере, когда распределение является центрированным и довольно узким).Таким образом, подпрограммы подбора в Scipy не могут определить, какие из них нужно изменить на любом данном шаге.
Исправление
Самый простой способ решить проблему - заблокировать один из ваших параметров.Вам не нужно менять уравнение, но вам нужно сделать хотя бы одну из a, x0 или sigma константой.Наилучший выбор параметра для исправления - это, вероятно, x0, поскольку определить среднее значение / медиану / режим ваших данных тривиально, просто получив координату x одного элемента данных, отличного от нуля в y.Вам также нужно будет немного более умным о том, как вы устанавливаете свои первоначальные догадки.Вот как это выглядит:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
x = np.linspace(0,1,20)
xdiff = x[1] - x[0]
y = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0])
# the mean/median/mode all occur at the x coordinate of the one datapoint that is non-zero in y
mean = x[np.argmax(y)]
# sigma should be tiny, since we want a narrow distribution
sigma = xdiff
# the scaling factor should be roughly equal to the "height" of the one datapoint
a = y.max()
def gaus(x,a,sigma):
return a*np.exp(-(x-mean)**2/(2*sigma**2))
bounds = ((1, .015), (20, 1))
popt,pcov = curve_fit(gaus, x, y, p0=[a, sigma], maxfev=20000, bounds=bounds)
residual = ((gaus(x,*popt) - y)**2).sum()
plt.figure(figsize=(8,6))
plt.plot(x,y,'b+:',label='data')
xdist = np.linspace(x.min(), x.max(), 1000)
plt.plot(xdist,gaus(xdist,*popt),'C0', label='fit distribution')
plt.plot(x,gaus(x,*popt),'ro:',label='fit')
plt.text(.1,6,"residual: %.6e" % residual)
plt.legend()
plt.show()
Вывод:
Лучшее исправление
Вы неТебе не нужно подходить, чтобы получить того гауссиана, которого ты хочешь.Вместо этого вы можете использовать простое выражение закрытой формы для вычисления необходимых параметров, как в функции fitonegauss в приведенном ниже коде:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def gauss(x, a, mean, sigma):
return a*np.exp(-(x - mean)**2/(2*sigma**2))
def fitonegauss(x, y, fwhm=None):
if fwhm is None:
# determine full width at half maximum from the spacing between the x points
fwhm = (x[1] - x[0])
# the mean/median/mode all occur at the x coordinate of the one datapoint that is non-zero in y
mean = x[np.argmax(y)]
# solve for sigma in terms of the desired full width at half maximum
sigma = fwhm/(2*np.sqrt(2*np.log(2)))
# max(pdf) == 1/(np.sqrt(2*np.pi)*sigma). Use that to determine a
a = y.max() #(np.sqrt(2*np.pi)*sigma)
return a, mean, sigma
N = 20
x = np.linspace(0,1,N)
y = np.zeros(N)
y[N//2] = 10
popt = fitonegauss(x, y)
plt.figure(figsize=(8,6))
plt.plot(x,y,'b+:',label='data')
xdist = np.linspace(x.min(), x.max(), 1000)
plt.plot(xdist,gauss(xdist,*popt),'C0', label='fit distribution')
residual = ((gauss(x,*popt) - y)**2).sum()
plt.plot(x, gauss(x,*popt),'ro:',label='fit')
plt.text(.1,6,"residual: %.6e" % residual)
plt.legend()
plt.show()
Вывод:
Преимуществ этого подхода много.Он намного эффективнее в вычислительном отношении, чем любая подгонка, он (по большей части) никогда не потерпит неудачу, и он дает вам гораздо больший контроль над фактической шириной дистрибутива, с которой вы в конечном итоге.
The *Функция 1041 * настроена так, что вы можете напрямую установить полную ширину на половину максимальной установленного распределения.Если вы оставите его неустановленным, код будет автоматически угадывать его по расстоянию между данными x.Похоже, это дает разумные результаты для вашего приложения.