Analyser des données expérimentales avec Sage
La densité optique est corrélée à la concentration de bactéries dans un liquide. Dans cet exemple, on réalise une régression linéaire sur une série de standards:
import numpy
var('OD, pente, ordonnee')
def droite_etalonnage(OD, pente, ordonnee):
"""Cette fonction décrit le modèle pour la droite d'étallonnage"""
return OD * pente + ordonnee
# Les données à utiliser sont sauvées dans des vecteurs Python
Concentration_cellules = numpy.array([2.60E+08, 3.14E+08, 3.70E+08, 4.62E+08, 8.56E+08, 1.39E+09, 1.84E+09])
densite_optique = numpy.array([0.083, 0.125, 0.213, 0.234,0.604, 1.092, 1.141])
OD_vs_concentration_cellules = zip(densite_optique, Concentration_cellules)
# On réalise la régression linéaire
params_regression = find_fit(OD_vs_concentration_cellules, droite_etalonnage, parameters=[pente, ordonnee], variables=[OD], initial_guess=[1e9, 3e8], solution_dict = True)
# On affiche les résultats de la régression linéaire
for param, valeur in params_regression.iteritems():
print('%s = %e' % (param, valeur))
# On dessine ensuite le graphique
data_plot = scatter_plot(OD_vs_Cell_concentration, markersize=20,facecolor='red', axes_labels=['OD a 600nm', 'Concentration de cellules (1/ml)'])
fit_plot = plot(droite_etalonnage(OD, params_regression[pente], params_regression[ordonnee]), (OD, 0, 1.2))
show(data_plot+fit_plot)
ordonnee = 1.237283e+08 pente = 1.324714e+09 ordonnee = 1.237283e+08 pente = 1.324714e+09 |
On peut ensuite visualiser la série de points expérimentaux:
# Les données mesurées sont les suivantes
temps_echantillons = numpy.array([0, 20, 40, 60, 80, 100, 120,140, 160, 180, 200, 220, 240, 280, 360, 380, 400, 420,440, 460, 500, 520, 540, 560,
580, 600, 620, 640, 660,680, 700, 720, 760, 1240, 1440, 1460, 1500, 1560])
densites_optiques = numpy.array([0.083, 0.087, 0.116, 0.119, 0.122,0.123, 0.125, 0.131, 0.138, 0.142, 0.158, 0.177, 0.213,0.234, 0.424, 0.604, 0.674, 0.726, 0.758, 0.828,
0.919,0.996, 1.024, 1.066, 1.092, 1.107, 1.113, 1.116, 1.12,1.129, 1.132, 1.135, 1.141, 1.109, 1.004, 0.984,0.972, 0.952])
concentrations = droite_etalonnage(densites_optiques, params_regression[pente], params_regression[ordonnee])
donnees_experimentales = zip(temps_echantillons, concentrations)
data_plot = scatter_plot(donnees_experimentales, markersize=20, facecolor='red', axes_labels=['temps (s)', 'Concentration de cellules (1/ml)'])
show(data_plot)
|
|
Maintenant, on recherche les paramètres d'un modèle de croissance sur ces données:
var('t, vitesse_max, temps_latence, y_max, y0')
def gompertz(t, vitesse_max, temps_latence, y_max, y0):
"""on définit un modèle de croissance selon Gompertz"""
return y0 + (y_max - y0) * numpy.exp(-numpy.exp(1.0 + vitesse_max * numpy.exp(1) * (temps_latence - t) / (y_max - y0)))
# On estime la valeur des paramètres à donner comme valeur initiale
vitesse_max_est = (1.4e9 - 5e8)/200.0
temps_latence_est = 100
y_max_est = 1.7e9
y0_est = 2e8
# On réalise ici une régression non linéaire
gompertz_params = find_fit(donnees_experimentales, gompertz, parameters=[vitesse_max, temps_latence, y_max, y0], variables=[t], initial_guess=[vitesse_max_est, temps_latence_est,
y_max_est, y0_est], solution_dict = True)
for param,value in gompertz_params.iteritems():
print(str(param) + ' = %e' % value)
y0 = 3.051120e+08 y_max = 1.563822e+09 temps_latence = 2.950639e+02 vitesse_max = 6.606311e+06 y0 = 3.051120e+08 y_max = 1.563822e+09 temps_latence = 2.950639e+02 vitesse_max = 6.606311e+06 |
On peut également comparer graphiquement le modèle avec les résultats expérimentaux
gompertz_model_plot = plot(gompertz(t, gompertz_params[vitesse_max], gompertz_params[temps_latence], gompertz_params[y_max],
gompertz_params[y0]), (t, 0, sample_times.max()))
show(gompertz_model_plot + data_plot)
|
|
