Ajustement non linéaire avec barres d'erreur pondérées - Minimizer/scipy.curve_fit/model.fit

Question

Je travaille sur un code d'adaptation Python pour Michaelies-Menten, une équation non linéaire, capable d'inclure des barres d'erreur pondérées. Pour le moment, j'ai essayé d'utiliser Minimizer et model.fit de lmfit. Bien que Minimizer ne comprend pas les barres d'erreur pondérées et model.fit semble être moins statistique que Minimizer.
Existe-t-il un moyen d'inclure des barres d'erreur pondérées dans Minimizer?
Est-ce que scipy.optimize.curve_fit serait un meilleur moyen d'adapter ce code?
Ou y at-il un autre programme d'adaptation qui serait mieux?

Mon code est ci-dessous

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 

Espérons que vous pouvez me aider! Vive

Answer 1

Si je comprends bien, vous voulez adapter le modèle décrit dans my_model aux données y (x) (dans les tableaux y et x) et utiliser l'incertitude y, std, de pondérer l'ajustement - réduisant au minimum (mode données)/incertitude plutôt que seulement données - modèle.

Pour ce faire, avec lmfit.Model, vous voulez passer un poids de 1./std (vérification probablement divde par zéro), comme:

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

(il n'a pas été clair pour moi pourquoi était à la fois std et std_1.

pour ce faire, avec Minimize, ajoutez std à la fcn_args tuple (arguments à transmettre à votre fonction objectif), et changer la fonction objectif de remplacer

return model - y 

avec

return (model -y)/std 

Avec cela, vous devriez être prêt à aller.

FWIW, Model.fit utilise Minimizer, donc ce n'est pas vraiment "moins statistique", c'est juste un accent différent. En outre, il existe probablement des moyens plus efficaces pour charger les données (peut-être une variation de numpy.loadtxt) mais ce n'est pas la question principale ici.