7

Pythonでyの不確実性を持ついくつかのデータポイントを当てはめようとしています。データは Python で x、y、yerr としてラベル付けされます。そのデータを loglog スケールで線形近似する必要があります。適合結果が適切かどうかの参考として、Python の結果と Scidavis の結果を比較します。

でcurve_fitを試しました

def func(x, a, b):
    return np.exp(a* np.log(x)+np.log(b))

popt, pcov = curve_fit(func, x, y,sigma=yerr)

kmpfitと同様に

def funcL(p, x):
   a,b = p
   return ( np.exp(a*np.log(x)+np.log(b)) )

def residualsL(p, data):
   a,b=p
   x, y, errorfit = data
   return (y-funcL(p,x)) / errorfit

a0=1
b0=0.1
p0 = [a0,b0]
fitterL = kmpfit.Fitter(residuals=residualsL, data=(x,y,yerr))
fitterL.parinfo = [{}, {}]
fitterL.fit(params0=p0)

そして、不確かさのないデータの 1 つ (つまり、yerr=1 を設定) でデータを当てはめようとすると、すべてがうまく機能し、結果は scidavis のものと同じになります。しかし、yerr をデータ ファイルの不確実性に設定すると、不穏な結果が得られます。Pythonでは、つまりa = 0.86になり、scidavisではa = 0.14になります。エラーが重みとして含まれていることについて何かを読みました。フィットを正しく計算するために、何かを変更する必要がありますか? または、私は何を間違っていますか？

編集：これはデータファイルの例です（x、y、yerr）

3.942387e-02    1.987800e+00    5.513165e-01
6.623142e-02    7.126161e+00    1.425232e+00
9.348280e-02    1.238530e+01    1.536208e+00
1.353088e-01    1.090471e+01    7.829126e-01
2.028446e-01    1.023087e+01    3.839575e-01
3.058446e-01    8.403626e+00    1.756866e-01
4.584524e-01    7.345275e+00    8.442288e-02
6.879677e-01    6.128521e+00    3.847194e-02
1.032592e+00    5.359025e+00    1.837428e-02
1.549152e+00    5.380514e+00    1.007010e-02
2.323985e+00    6.404229e+00    6.534108e-03
3.355974e+00    9.489101e+00    6.342546e-03
4.384128e+00    1.497998e+01    2.273233e-02

そして結果：

in python: 
   without uncertainties: a=0.06216 +/- 0.00650 ; b=8.53594 +/- 1.13985
   with uncertainties: a=0.86051 +/- 0.01640 ; b=3.38081 +/- 0.22667 
in scidavis:
   without uncertainties: a  = 0.06216 +/- 0.08060; b  = 8.53594 +/- 1.06763
   with uncertainties: a  = 0.14154 +/- 0.005731; b  = 7.38213 +/- 2.13653
4

1 に答える 1

3

私は何かを誤解しているに違いない。投稿されたデータは次のようには見えません

f(x,a,b) = np.exp(a*np.log(x)+np.log(b))

赤い線は の結果scipy.optimize.curve_fit、緑の線は scidavis の結果です。

私の推測では、どちらのアルゴリズムも適切な適合に向かって収束していないため、結果が一致しないことは驚くべきことではありません。

scidavis がそのパラメーターを見つける方法を説明することはできませんが、私が理解している定義によれば、scipyよりも最小二乗残差が小さいパラメーターを見つけていscidavisます。

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as optimize

def func(x, a, b):
    return np.exp(a* np.log(x)+np.log(b))

def sum_square(residuals):
    return (residuals**2).sum()

def residuals(p, x, y, sigma):
    return 1.0/sigma*(y - func(x, *p))

data = np.loadtxt('test.dat').reshape((-1,3))
x, y, yerr = np.rollaxis(data, axis = 1)
sigma = yerr

popt, pcov = optimize.curve_fit(func, x, y, sigma = sigma, maxfev = 10000)
print('popt: {p}'.format(p = popt))
scidavis = (0.14154, 7.38213)
print('scidavis: {p}'.format(p = scidavis))

print('''\
sum of squares for scipy:    {sp}
sum of squares for scidavis: {d}
'''.format(
          sp = sum_square(residuals(popt, x = x, y = y, sigma = sigma)),
          d = sum_square(residuals(scidavis, x = x, y = y, sigma = sigma))
      ))

plt.plot(x, y, 'bo', x, func(x,*popt), 'r-', x, func(x, *scidavis), 'g-')
plt.errorbar(x, y, yerr)
plt.show()

収量

popt: [ 0.86051258  3.38081125]
scidavis: (0.14154, 7.38213)
sum of squares for scipy:    53249.9915654
sum of squares for scidavis: 239654.84276

ここに画像の説明を入力

于 2012-08-06T18:20:54.293 に答える