我试图用高斯函数拟合一些非常清晰的数据,但由于某种原因,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 ]
而剧情:
编辑: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 /错误:我发现这两个数组分别是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')
工作正常...
1条答案
按热度按时间gc0ot86w1#
答:为发现此问题的任何人提供文档,显然这(我认为)是scipy中的一个bug。我发现两个数组是np.float32和np.float64,显然scipy无法正确处理这一问题,如果我这样做:
则装配失败
工作,本质上的问题是最小二乘最小化不能处理不同精度的输入。12显示器