2016-04-05 3 views
0

データポイントに合っています。エラーのないフィッティングが楽観的ではないように見えるので、今度は、エラーを実装するデータを各ポイントに合わせようとしています。マイフィット機能は以下の通りです:エラーバーを使用してScipyカーブフィッティングを実行し、フィッティングパラメータの標準エラーを取得する方法?

def fit_func(x,a,b,c): 
    return np.log10(a*x**b + c) 

その後、私のデータポイントは以下の通りです:

r = [ 0.00528039,0.00721161,0.00873037,0.01108928,0.01413011,0.01790143,0.02263833, 0.02886089,0.03663713,0.04659512,0.05921978,0.07540126,0.09593949, 0.12190075,0.15501736,0.19713563,0.25041524,0.3185025,0.40514023,0.51507869, 0.65489938,0.83278859,1.05865016,1.34624082] 
logf = [-1.1020581079659384, -1.3966927245616112, -1.4571368537041418, -1.5032694247562564, -1.8534775558300272, -2.2715812166948304, -2.2627690390113862, -2.5275290780299331, -3.3798813619309365, -6.0, -2.6270989211307034, -2.6549656159564918, -2.9366845162570079, -3.0955026428779604, -3.2649261507250289, -3.2837123017838366, -3.0493752067042856, -3.3133647996463229, -3.0865051494299243, -3.1347499415910169, -3.1433062918466632, -3.1747394718538979, -3.1797597345585245, -3.1913094832146616] 

私のデータは、対数スケール、logfにあるので、その後、各データポイントのエラーバーは左右対称ではありません。上部のエラーバーと下のエラーバーは以下の通りです:

upper = [0.070648916083227764, 0.44346256268274886, 0.11928131794776076, 0.094260899008089094, 0.14357124858039971, 0.27236750587684311, 0.18877122991380402, 0.28707938182603066, 0.72011863806906318, 0, 0.16813325716948757, 0.13624929595316049, 0.21847915642008875, 0.25456116079315372, 0.31078368240910148, 0.23178227464741452, 0.09158189214515966, 0.14020538489677881, 0.059482730164901909, 0.051786777740678414, 0.041126467609954531, 0.034394612910981337, 0.027206248503368613, 0.021847333685597548] 
lower = [0.06074797748043137, 0.21479225959441428, 0.093479845697059583, 0.077406149968278104, 0.1077175009766278, 0.16610073183912188, 0.13114254113054535, 0.17133966123838595, 0.57498950902908286, 2.9786837094190934, 0.12090437578535695, 0.10355760401838676, 0.14467588244034646, 0.15942693835964539, 0.17929440903034921, 0.15031667827534712, 0.075592499975030591, 0.10581886912443572, 0.05230849287772843, 0.04626422871423852, 0.03756658820680725, 0.03186944137872727, 0.025601929615431285, 0.02080073540367966] 

私のようにフィッティングを持っている:

popt, pcov = optimize.curve_fit(fit_func, r, logf,sigma=[lower,upper]) 
logf_fit = fit_func(r,*popt) 

しかし、これは間違っている、どのように私は、上側が含まれるようにscipyのダウンロードからカーブフィッティングを実装することができますし、より低いエラー?どのようにフィッティングパラメータa、b、cのフィッティング誤差を得ることができますか?

+0

なぜあなたは '' a * x ** b + c''にすぐにフィットし、エラーを対称に保つのですか? – tBuLi

+0

ありがとうございます。私はこの機能を試しましたが、それは私に正しい物理学を与えません。それから、私は上記のアプローチを試しました。物理的に正しいものになります。 –

+0

[pythonのoptimize.leastsqメソッドを使って適合するパラメータの標準エラーを取得する](http://stackoverflow.com/questions/14581358/getting-standard-errors-on-fitted-parameters-using-the-optimize -leastsq-method-i) –

答えて

1

カスタム重みでscipy.optimize.leastsqを使用することができます。

import scipy.optimize as optimize 
import numpy as np 
# redefine lists as array 
x=np.array(r) 
y=np.array(logf) 
errup=np.array(upper) 
errlow=np.array(lower) 
# error function 
def fit_func(x,a,b,c): 
    return np.log10(a*x**b + c) 
def my_error(V): 
    a,b,c=V 
    yfit=fit_func(x,a,b,c) 
    weight=np.ones_like(yfit) 
    weight[yfit>y]=errup[yfit>y] # if the fit point is above the measure, use upper weight 
    weight[yfit<=y]=errlow[yfit<=y] # else use lower weight 
    return (yfit-y)**2/weight**2 
answer=optimize.leastsq(my_error,x0=[0.0001,-1,0.0006]) 
a,b,c=answer[0] 
print(a,b,c) 

それは動作しますが、間違ったドメイン(負の数)に行くことができるログがあるので、初期値に非常に敏感であるし、それが失敗しました。ここで私はかなり近いデータであるa=9.14464745425e-06 b=-1.75179880756 c=0.00066720486385を見つける。

+0

ありがとうございます。フィッティングパラメータa、b、cの不確実性をどのように出力できるか知っていますか? –

+0

また、my_error関数の戻り文が(abs(yfit-y)/ weight)** 2でなければならないのでしょうか?それは最小二乗であるので、体重の四角だけではありません。 –

+0

はい、そうです、正方形でコードを修正しました。結果はより良いです。不確実性については、[doc](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html)を参照してください。特に、より多くのオプションを返す 'full_output'オプションがありますあなたが不確実性を推定することができる 'cov_x'は、私は(このドメインではあまりよくありません)と思います。 – JPG

関連する問題