2

I am working on a Python fitting code for Michaelies-Menten, a non-linear equation, able to include weighted errorbars. At the moment I have tried using Minimizer and model.fit from lmfit. Though Minimizer does not include weighted errorbars and model.fit seems to be less statistical than Minimizer.
Is there a way to include weighted errorbars in Minimizer?
Would scipy.optimize.curve_fit be a better way to fit this code?
Or is there another fitting program that would be better?

My code is below

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

Hope you can help me! Cheers

4

1 回答 1

0

如果我理解正确，您希望将 y(x) 中描述的模型拟合到数据y ( x )my_model（在数组y和- 模型。xystd

要做到这一点lmfit.Model，你需要传递一个权重1./std（可能检查 divde-by-zero），如下所示：

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

（我不清楚为什么同时存在std和std_1。

为此，请Minimize添加std到fcn_args元组（要传递给目标函数的参数），并更改目标函数以替换

return model - y

和

return (model -y)/std

有了这个，你应该准备好了。

FWIW, Model.fituses Minimizer，所以它并不是真正的“较少统计”，它只是一个不同的重点。

顺便说一句，可能有更有效的方法来加载数据（可能是 的一些变体numpy.loadtxt），但这不是这里的主要问题。

于 2017-05-11T21:45:11.663 回答