Weighted Errorbarsを使用した非線形フィッティング - Minimizer/scipy.curve_fit/model.fit

Question

Michaelies-MentenのPythonフィッティングコード、非線形方程式に加重誤差バーを含めることができます。現時点では、私はMinimizerとmodel.fitをlmfitから使ってみました。 Minimizerには加重誤差バーが含まれていませんが、model.fitはMinimizerよりも統計的には小さいようです。
重み付きエラーバーをMinimizerに含める方法はありますか？
scipy.optimize.curve_fitはこのコードに適した方法でしょうか？
もう1つのフィッティングプログラムがありますか？あなたが私を助けることを願って

def michealies_menten(path, conc, substrate): 
os.chdir(path) 
h = open('average_STD_data.dat', 'r') 
f = open('concentration.dat', 'r') 

x = [] 
y = [] 
std = [] 
std_1 = [] 

for line in h: 
    line = line.split(',') 
    y.append(line[0]) 
    std.append(line[1]) 

for line in f: 
    x = line.split(',') 

for i in range(len(std)): 
    new = 1.0/(float(std[i])**2.0) 
    std_1.append(float(new)) 

std.insert(0, '0') 
std_1.insert(0, '0') 
x.insert(0, '0') 
y.insert(0, '0') 
y=map(float,y) 
x=map(float,x) 
std=map(float,std) 
std_1=map(float,std_1) 

x=np.array(x) 
y=np.array(y) 
std_1=np.array(std_1) 

####Model.fit code: 
def my_model(x, Vmax, Km): 
    return Vmax * x/(Km + x) 
gmodel = Model(my_model) 
result = gmodel.fit(y, x=x, Vmax=4000.0, Km=3.0, weights=std_1) 
print result.fit_report() 
Vmax_1=result.params['Vmax'].value 
Km_1=result.params['Km'].value 

model = (Vmax_1*x/(Km_1+x)) 



###Minimizer code: 
def get_residual(params, x, y): 
    Vmax = params['Vmax'] 
    Km = params['Km'] 
    model = Vmax * x/(Km + x) 
    return model - y 

#Parameters definition for LmFit 
params = Parameters() 
params.add('Vmax',value=4000., min=0) 
params.add('Km',value=3., min=0) 

#Produces the Km and Vmax value which is then extranted 
minner = Minimizer(get_residual, params, fcn_args=(x,y)) 
result = minner.minimize() 
print "Result of minization, deviation from y value:", result.residual 

#Resulting in the final y-data, which gives the fitted data. 
final = y + result.residual 
print "I am here, Final:", final 

#Gives report over the minimize function 
print "Result.params:" 
result.params.pretty_print() 
print "Report_fit:" 
report_fit(result) 

#Transfer lmFit output of the minimize(result) to variable for further use 
Vmax=result.params['Vmax'].value 
Km=result.params['Km'].value 
print "Fitted - Heres Km", Km 
print "Fitted - Heres Vmax", Vmax 

#Draw the different graphs 
#plt.plot(x,fitfunc_michment(x, *popt),'b',label='curve_fit') 
plt.plot(x,final,'r',label='lmfit') 
plt.plot(x,model,'b',label='model') 
plt.plot(x,y,'rs',label='Raw') 
plt.errorbar(x,y,yerr=std, ecolor='r', linestyle="None") 
plt.xlabel(s.join(['Concentration of ', subst ,' [', conc_unit,']']),fontsize=12) 
plt.ylabel('Intensity [A.U.]', fontsize=12) 
plt.savefig('Michaelis_Menten_plot.png', bbox_inches='tight') 
plt.show() 
print 'Hello World, i am the Km value: ', Km 
print 'Vmax value: ', Vmax 

の下

私のコードです！ 乾杯

Answer 1

私が正しく理解していれば、あなたは（X）（y配列にし、x）データyにmy_modelで説明したモデルに適合し、フィット感を重み付けする、y、stdに不確実性を使用したい - 最小化（データモード）/不確定性ではなく、データモデルです。 lmfit.Modelでこれを行うには

は、あなたが1./stdの重量に渡したい（おそらくdivdeゼロによるのチェック）、のように：

result = gmodel.fit(y, x=x, Vmax=4000.0, Km=3.0, weights=1.0/(std)) 

（それはなぜそこに私には明確ではありませんでしたstdとstd_1でもあった。

を置き換えるために、 Minimizeでこれを行う（引数は、あなたの目的関数に渡される） fcn_argsタプルのよう stdを追加し、目的関数を変更するにはそれと

return (model -y)/std 

と

return model - y 

、あなたは行く準備ができなければなりません。

FWIW、Model.fitはMinimizerを使用しているため、統計的にはそれほど重要ではありません。

これ以外にも、おそらくデータを読み込むための効率的な方法があります（おそらくnumpy.loadtxtのバリエーション）が、これは主な質問ではありません。