3  Ajustament logístic

3.1 Ajustament de dades

En aquesta pràctica anem a simular els exemples de Gause a la dècada dels 30 que va intentar reproduir experimentalment el creixement logístic d’algunes espècies. Concretament seguirem l’experiment sobre un fong del llevat, Saccharomyces cerevisiae. El Saccharomyces cerevisiae és un organisme que forma part del grup de llevats de gemmació. Probablement es tracta del tipus de llevat més important de la història, ja que des dels temps més remots s’ha emprat tant en fermentacions de begudes i aliments com la cervesa o el pa. Es creu que aquest llevat fou originàriament aïllat de la pell de raïm (hom troba aquest llevat com a component de la fina capa blanca a la pell d’algunes fruites fosques com les prunes; es troba entre la cera de la cutícula). Mantingut a unes condicions constants, Gause va obtenir les següents dades experimentals:

Temps	Volum
0	0.37
1.5	1.63
9	6.2
10	8.87
18	10.66
18	10.97
23	12.5
25.5	12.6
27	12.9
34	13.27
38	12.77
42	12.87
45.5	12.9
47	12.7

Ho escrivim en un fitxer “cervisiae.txt” i ho llegim

from numpy import *
from scipy.io import write_array
from scipy.io import read_array
#Les llegim:
data = read_array(file("cervisiae.txt"), separator=" ", columns=(0,1))

Quan ho tenim ho representem:

from pylab import *
x=data[:,0]
y=data[:,1]
plot(x,y, 'o--')
show()

i obtenim una gràfica semblant. Anem a ajustar-hi una logística. En primer lloc recordem que l’equació logística, en aquestes variables, llegeix: \[y'(x)= r y \left(1- \frac{y}{K} \right),\] que té com a solució: \[y(x)= \frac{y(0)K}{y(0)+(K-y(0))e^{-rx}}\] Anem a calcular-ho per mitjà dels mínims quadrats. Per això usarem la instrucció leastsq que ens permet optimitzar una funció amb paràmetres:

from scipy import *
# Fit the first set
fitfunc = lambda p, x: p[0]*p[1]/(p[0]+(p[1]-p[0])*exp(-p[2]*x))  # funcio a fitar
errfunc = lambda p, x, y: fitfunc(p,x) -y          # Distancia a la funcio a fitar
p0 = [0.7,12.9,0.25]                        # Prova inicial.
p1,success = optimize.leastsq(errfunc, p0[:], args = (x, y))
temps=linspace(x[0],x[-1],100)
plot(x,y,"ro")+plot(temps,fitfunc(p1,temps),"--")        # Ho dibuixem tot

La qüestió és, com podem trobar unes aproximacions per als paràmetres \(y(0)\), \(r\) i \(K\)? \[y(0) = 0.37 M \simeq 12.7 r \simeq y'(0)/y(0)\] i així podem fer:

derivada1= (y[1]-y[0])/(x[1]-x[0])
p0 = [x[0],y[-1],derivada1/y[0]]
p1,success = optimize.leastsq(errfunc, p0, args = (x, y))

Proveu de fer el mateix amb les dades de Carlson (1913) en estudiar un altre tipus de llevat:

Temps	Volum	Temps	Volum	Temps	Volum
1	9.6	7	174.6	13	594.8
2	18.3	8	257.3	14	629.4
3	29.0	9	350.7	15	640.8
4	47.2	10	441.0	16	651.1
5	71.1	11	513.3	17	655.9
6	119.1	12	559.7	18	659.6

Feu el mateix amb un ajustament Gomperzià: \[y'(x)= r y \log\left(K/y(x)\right)\] que té com a solució: \[y(x)= \left(\frac{y(0)}{K}\right)^{e^{-rt}}\]