python - scipy.optimize.leastsq 返回最佳猜测参数而不是新的最适合

Question

我想将洛伦兹峰拟合到一组数据 x 和 y，数据很好。 OriginLab 等其他程序非常适合它，但我想使用 python 自动进行拟合，所以我有以下基于 http://mesa.ac.nz/?page_id=1800 的代码

我遇到的问题是 scipy.optimize.leastsq 返回我传递给它的相同初始猜测参数的最佳拟合，实际上什么都不做。这是代码。

#x, y are the arrays with the x,y  axes respectively
#defining funcitons
def lorentzian(x,p):
  return p[2]*(p[0]**2)/(( x - (p[1]) )**2 + p[0]**2)

def residuals(p,y,x):
  err = y - lorentzian(x,p)
  return err      

p = [0.055, wv[midIdx], y[midIdx-minIdx]]   
pbest = leastsq(residuals, p, args=(y, x), full_output=1)
best_parameters = pbest[0]
print p
print pbest

p 是初始猜测，best_parameters 是从 leastsq 返回的“最适合”参数，但它们始终相同。

这是 full_output=1 返回的结果(长数值数组已经缩短但仍然具有代表性)

    [0.055, 855.50732, 1327.0]
    (array([  5.50000000e-02,   8.55507324e+02,   1.32700000e+03]), 
    None, {'qtf':array([ 62.05192947,  69.98033905,  57.90628052]), 
    'nfev': 4, 
    'fjac': array([[-0.,  0.,  0., 0.,  0.,  0.,  0.,],
    [ 0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
    0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
    0.,  0.],
    [ 0.,  0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
    0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
    0.,  0.]]), 
    'fvec': array([  62.05192947,   69.98033905,   
    53.41218567,   45.49879837,   49.58242035,   36.66483688,
    34.74443436,   50.82238007,   34.89669037]), 
    'ipvt': array([1, 2, 3])},  
    'The cosine of the angle between func(x) and any column of the\n  Jacobian 
    is at most 0.000000 in absolute value', 4)

谁能看出哪里出了问题？

Answer 1

快速谷歌搜索提示数据是单精度的问题(您的其他程序几乎肯定也升级为 double ，尽管这显然也是 scipy 的问题，另请参见 bug report )。如果您查看 full_output=1 结果，您会发现雅可比矩阵在任何地方都近似为零。因此，明确给出 Jacobian 行列式可能会有所帮助(尽管那样你可能想要向上转换，因为你可以用单精度获得的相对误差的最小精度非常有限)。

解决方案:最简单且在数值上最好的解决方案(当然给出真正的雅可比行列式也是一种奖励)是只转换你的 x 和 y 数据到 double (例如 x = x.astype(np.float64) 会做)。

我不建议这样做，但您也可以通过手动设置 epsfcn 关键字参数(可能还有 tolerance 关键字参数)来修复它，epsfcn=np .finfo(np.float32).eps。这似乎在某种程度上解决了这个问题，但是(因为大多数计算都是用标量计算的，并且标量不会强制在你的计算中进行向上转换)计算是在 float32 中完成的，精度损失似乎相当大，至少当不是提供 Dfunc。