Я пытаюсь подогнать гауссиан к каким-то довольно четким данным, и по какой-то причине scipy.optimize.curve_fit совершенно не меняет параметры по сравнению с первоначальным предположением, я не знаю, почему.Вот фрагмент кода:
def gauss(x, A, mu, sigma):
return A/(sigma*np.sqrt(2*np.pi))*np.exp(-((x-mu)**2/(2*sigma**2)))
p0 = [grad[ir_max]/5, plot_r[ir_max], 0.05]
print(p0)
fit, pcov = curve_fit(gauss, plot_r[ir_max-dr:ir_max+dr], grad[ir_max-
dr:ir_max+dr], p0=p0)
print(fit)
print(plot_r[ir_max-dr:ir_max+dr])
print(grad[ir_max-dr:ir_max+dr])
plot(plot_r[ir_max-dr:ir_max+dr], grad[ir_max-dr:ir_max+dr], 'bo')
plot_r2 = np.linspace(plot_r[ir_max-dr], plot_r[ir_max+dr], 100)
plot(plot_r2, gauss(plot_r2, *fit), 'r--')
Индексы должны взять подмножество большего набора данных, вот вывод для довольно ограниченной выборки вокруг пика:
[7.7160651775860245, 1.641777, 0.05]
[7.71606518 1.64177704 0.05 ]
[1.5620524 1.5779973 1.5939423 1.6098871 1.6258321 1.641777 1.657722
1.673667 1.6896119 1.7055569]
[ 7.21488949 15.13438187 25.0198808 33.35524257 37.91767649 38.58032589
35.52668657 28.27396106 18.12926291 9.1928141 ]
Иplot:
Редактировать: scipy.omtimize.leastsq также имеет аналогичную проблему, но выполнение минимизации наименьших квадратов вручную с помощью scipy.optimize.minimize, кажется, работаетОК:
def func(p):
A, mu, sigma = p
resid = gauss(plot_r[ir_max-dr:ir_max+dr], A, mu, sigma)
return grad[ir_max-dr:ir_max+dr] - resid
from scipy.optimize import leastsq
fit2 = leastsq(func, p0, full_output=True)
print(fit2[0])
print(fit2)
def func2(p):
A, mu, sigma = p
resid = gauss(plot_r[ir_max-dr:ir_max+dr], A, mu, sigma)
return np.sum( np.power(grad[ir_max-dr:ir_max+dr] - resid,2) )
from scipy.optimize import minimize
fit3 = minimize(func2, p0, method='Nelder-Mead')
print(fit3.x)
plot(plot_r[ir_max-dr:ir_max+dr], grad[ir_max-dr:ir_max+dr], 'bo')
plot_r2 = np.linspace(plot_r[ir_max-dr], plot_r[ir_max+dr], 100)
plot(plot_r2, gauss(plot_r2, *fit), 'r--')
plot(plot_r2, gauss(plot_r2, *fit2[0]), 'm--')
plot(plot_r2, gauss(plot_r2, *fit3.x), 'c--')
с результатами print операторов:
[7.71606518 1.64177704 0.05 ]
(array([7.71606518, 1.64177704, 0.05 ]), None, {'fvec': array([-10.05377583, -12.15582504, -13.93739459, -16.87939946,
-20.59538122, -22.98496647, -22.98637114, -21.96068096,
-20.82792474, -18.09731079]), 'nfev': 21, 'fjac': array([[ 9.57864291e+03, 0.00000000e+00, 0.00000000e+00,
5.34522484e-01, -0.00000000e+00, -0.00000000e+00,
-0.00000000e+00, 5.34522484e-01, -0.00000000e+00,
2.67261242e-01],
[ 0.00000000e+00, -0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, -0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00]]), 'ipvt': array([3, 2, 1], dtype=int32), 'qtf': array([ 34.78215243, 6.75567556, -13.93739459])}, 'The relative error between two consecutive iterates is at most 0.000000', 2)
[4.21310821 1.6364875 0.04207018]
и сюжетом: 
РЕДАКТИРОВАТЬ2 / BUG: Документирование для любого, кто нашел это, по-видимому это (то, что я думаю) ошибка в scipy.Я обнаружил, что эти два массива были np.float32 и np.float64, и, очевидно, scipy не может справиться с этим должным образом, если я это сделаю:
xdata = np.copy(plot_r[ir_max-dr:ir_max+dr])
ydata = np.copy(grad[ir_max-dr:ir_max+dr])
fit, pcov = curve_fit(gauss, xdata, ydata, p0=p0, method='lm')
подгонка не удалась, но
xdata = xdata.astype(dtype=np.float32)
ydata = ydata.astype(dtype=np.float32)
fit, pcov = curve_fit(gauss, xdata, ydata, p0=p0, method='lm')
работает нормально...