1

我一直在尝试用 scipy.optimize.curve_fit 拟合一些直方图数据，但到目前为止，我还没有一次能够生成与我的猜测参数显着不同的拟合参数。

如果我发现拟合中的更多神秘参数陷入局部最小值，我不会感到非常惊讶，但即使是线性系数也不会偏离我最初的猜测！

如果你以前见过这样的事情，我会喜欢一些建议。最小二乘最小化例程是否不适用于某些类型的函数？

我试试这个

import numpy as np
from matplotlib.pyplot import *
from scipy.optimize import curve_fit

def grating_hist(x,frac,xmax,x0):
  #  model data to be turned into a histogram
  dx = x[1]-x[0]
  z = np.linspace(0,1,20000,endpoint=True)
  grating = np.cos(frac*np.pi*z)
  norm_grating = xmax*(grating-grating[-1])/(1-grating[-1])+x0
  # produce the histogram
  bin_edges = np.append(x,x[-1]+x[1]-x[0])
  hist,bin_edges = np.histogram(norm_grating,bins=bin_edges)
  return hist

x = np.linspace(0,5,512)
p_data = [0.7,1.1,0.8]
pct = grating_hist(x,*p_data)
p_guess = [1,1,1]
p_fit,pcov = curve_fit(grating_hist,x,pct,p0=p_guess)

plot(x,pct,label='Data')
plot(x,grating_hist(x,*p_fit),label='Fit')
legend()

show()

print 'Data Parameters:', p_data
print 'Guess Parameters:', p_guess
print 'Fit Parameters:', p_fit
print 'Covariance:',pcov

我看到了这个：http: //i.stack.imgur.com/GwXzJ.png（我是新来的，所以我不能发布图片）

Data Parameters: [0.7, 1.1, 0.8]
Guess Parameters: [1, 1, 1]
Fit Parameters: [ 0.97600854  0.99458336  1.00366634]
Covariance: [[  3.50047574e-06  -5.34574971e-07   2.99306123e-07]
 [ -5.34574971e-07   9.78688795e-07  -6.94780671e-07]
 [  2.99306123e-07  -6.94780671e-07   7.17068753e-07]]

哇？我很确定这不是 xmax 和 x0 变化的局部最小值，而且距离全局最小值最佳拟合还有很长的路要走。即使有更好的猜测，拟合参数仍然不会改变。曲线函数的不同选择（例如两个正态分布的总和）确实会为相同的数据产生新的参数，所以我知道这不是数据本身。为了以防万一，我也对 scipy.optimize.leastsq 本身进行了同样的尝试，但没有骰子；参数仍然不动。如果您对此有任何想法，我很想听听他们的意见！

4

2 回答 2

1

您面临的问题实际上不是由于curve_fit（或leastsq）。这是由于您的优化问题目标的情况。在您的情况下，目标是您试图最小化的残差平方的总和。现在，如果您在初始条件附近仔细观察您的目标，例如使用下面的代码，它只关注第一个参数：

p_ind = 0
eps = 1e-6
n_points = 100

frac_surroundings = np.linspace(p_guess[p_ind] - eps, p_guess[p_ind] + eps, n_points)

obj = []
temp_guess = p_guess.copy()
for p in frac_surroundings:
    temp_guess[0] = p
    obj.append(((grating_hist(x, *p_data) - grating_hist(x, *temp_guess))**2.0).sum())

py.plot(frac_surroundings, obj)
py.show()

您会注意到景观是分段常数（您可以轻松检查其他参数的情况是否相同。问题在于这些部分的数量级为 10^-6，而拟合过程的初始步骤大约在 10^-8 左右，因此该过程很快结束，得出的结论是您无法从给定的初始条件改进。您可以尝试通过更改epsfcn参数来修复它curve_fit，但您会很快注意到景观，除了分段常数, 也非常“坚固”。换句话说，curve_fit它根本不适合这种问题，这对于基于梯度的方法来说简直是困难的，因为它是高度非凸的。也许，一些随机优化方法可以做得更好. 然而，这是一个不同的问题/问题。

于 2021-03-04T01:54:53.930 回答
0

我认为这是一个局部最小值，或者算法因非平凡原因而失败。将数据拟合到输入要容易得多，而不是将数据的统计描述拟合到输入的统计描述。

这是这样做的代码的修改版本：

z = np.linspace(0,1,20000,endpoint=True)

def grating_hist_indicator(x,frac,xmax,x0):
  #  model data to be turned into a histogram
  dx = x[1]-x[0]
  grating = np.cos(frac*np.pi*z)
  norm_grating = xmax*(grating-grating[-1])/(1-grating[-1])+x0
  return norm_grating

x = np.linspace(0,5,512)
p_data = [0.7,1.1,0.8]
pct = grating_hist(x,*p_data)

pct_indicator = grating_hist_indicator(x,*p_data)
p_guess = [1,1,1]
p_fit,pcov = curve_fit(grating_hist_indicator,x,pct_indicator,p0=p_guess)

plot(x,pct,label='Data')
plot(x,grating_hist(x,*p_fit),label='Fit')
legend()
show()
于 2013-06-28T11:52:59.220 回答