4

I want to fit a lorentzian peak to a set of data x and y, the data is fine. Other programs like OriginLab fit it perfectly, but I wanted to automate the fitting with python so I have the below code which is based on http://mesa.ac.nz/?page_id=1800

The problem I have is that the scipy.optimize.leastsq returns as the best fit the same initial guess parameters I passed to it, essentially doing nothing. Here is the code.

#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 are the initial guesses and best_parameters are the returned 'best fit' parameters from leastsq, but they are always the same.

this is what returned by the full_output=1 (the long numeric arrays have been shortened but are still representitive)

    [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)

can anyone see whats wrong?

4

1 に答える 1

5

簡単なGoogle検索は、データが単精度であるという問題を示唆しています（他のプログラムもほぼ確実に倍精度にアップキャストされますが、これは明らかにscipyの問題でもあります。このバグレポートも参照してください）。結果を見るfull_output=1と、ヤコビアンがどこでもゼロに近似されていることがわかります。そのため、ヤコビアンを明示的に与えると役立つ場合があります (ただし、単精度で取得できる相対誤差の最小精度は非常に限られているため、アップキャストすることもできます)。

解決策：最も簡単で数値的に最良の解決策（もちろん、実際のヤコビアンを与えることもおまけです）は、データxとyデータを倍精度にキャストするx = x.astype(np.float64)ことです（たとえば、実行します）。

epsfcn私はこれをお勧めしませんが、キーワード引数 (およびおそらく許容キーワード引数) を手動で設定することで修正できる場合がありますepsfcn=np.finfo(np.float32).eps。これはある意味で問題を解決しているように見えますが、(ほとんどの計算はスカラーを使用しており、スカラーは計算でアップキャストを強制しないため) 計算は float32 で行われ、少なくともそうでない場合は精度の損失がかなり大きいようです。 Dfunc を提供します。

于 2012-09-18T12:36:27.487 に答える