scipy.optimize.leastsq返回最佳猜测参数,而不是新的最佳拟合

Question

我想将lorentzian峰值拟合到一组数据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

一个快速的谷歌搜索暗示数据是单精度的问题(你的其他程序几乎肯定是向上翻倍到双精度,虽然这显然也是scipy的问题,另见这个错误报告).如果你看一下你的full_output=1结果,你就会看到Jacobian到处都是零.因此,明确地给予雅可比行星可能有所帮助(尽管即使这样你也可能想要向上转换,因为单精度可以得到的相对误差的最小精度非常有限).

解决方案:最简单和数值最佳的解决方案(当然,给予真正的雅可比行列式也是一种奖励)就是将您的数据x和y数据转换为双精度(x = x.astype(np.float64)例如).

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