Fitting e modelli matematici


Calibrazione per un modello di crescita di parameci aurelia

Partiamo dai seguenti dati per una popolazione di parameci aurelia che inseriamo in un foglio elettronico di Excel
(che si puo' scaricare qui)


A
B
C
D
E
F
G
1
giorni
popolazione





2
2
25





3
3
51





4
4
67





5
5
98





6
6
197





7
7
287





8
8
339





9
9
409





10
10
512





11
11
519





12
12
534





13
13
544





14
14
551





15
15
556






Come prima operazione, vogliamo effettuare un fitting con il modello logistico di crescita mediante trasformazione delle variabili (in modo a ridursi a una equazione lineare) e poi mediante regressione lineare.

Alcuni modelli e trasformazioni linearizzanti:

Esponenziale: y=a exp(bt) ------> log(y)=log(a)+bt
Logistica:
y=K/(1+a exp(-bt)) ------> log(K/y-1)=log(a)-bt
Sigmoide:
y=K/(1+a tb)
------> log(K/y-1)=log(a)+b*log(t)


La funzione logistica e' y=K/(1+a exp(-bt)). Dai dati supponiamo che la popolazione limite (di equilibrio) sia K=560.  Abbiamo allora  K/y-1= a exp(-bt), e presi i logaritmi di ambo i membri, log(K/y-1)=log(a)-bt. Determiniamo allora a e b mediante regressione lineare. Inseriamo nella terza colonna i dati calcolando log(K/y-1).  A tale scopo scriviamo della cella C2 la formula =LOG((560/$B$2-1);2,718281) etc.
Otteniamo:


A
B
C
D
E
F
G
1
giorni
popolazione
log(K/y-1)




2
2
25
3,063391856




3
3
51
2,300623085



4
4
67
1,995817163



5
5
98
1,550597885



6
6
197
0,611199292



7
7
287
-0,050010436



8
8
339
-0,427837536




9
9
409
-0,996435623




10
10
512
-2,367124336




11
11
519
-2,53833259




12
12
534
-3,022300222




13
13
544
-3,526361599




14
14
551
-4,114511486




15
15
556
-4,934475437






Disegnamo ora il grafico con l'opzione "chart" inserendo come asse x i giorni e asse y i dati della colonna C. Al tempo stesso inseriamo la retta di regressione ("trendline" - "linea di tendenza" nella versione italiana) insieme con la sua equazione.





Dunque  log(a)=4,332, a=76,0963, b=0,6143 Si trova la curva logistica y= 560*1/(1+76.0963*exp(-.6134*t)).

In alternativa si puo' effettuare mediante lo strumento "analisi di dati". (nella versione inglese si trova sotto il menu "Tools" e si chiama "Data analysis".)

Nota: Potrebbe succedere che non si trovi tale opzione tra gli strumenti. Bisogna in primo luogo cercare tra gli "Add-ins" ("Aggiunte") e -se lo si trova- aggiungere "analisi di dati". Potrebbe pero' non trovarsi. In questo caso tale aggiunta non e' stata installata (con l'installazione standard di "Office" non viene installata): pertanto e' necessario usare il CD di installazione.


Vediamo ora come possiamo ottimizzare la scelta dei parametri usando il "Solver" ("Risolutore") di Excel. (Come nel caso di "analisi di dati", p
otrebbe succedere che non si trovi tale opzione tra gli strumenti. Si proceda come nella nota precedente.)

Nella seconda colonna inseriamo i dati in base alla formula data dal modello di crescita logistica  y=K/(1+a exp(-bt)). Per fare questo, si inserisce nella terza colonna (C) la formula scrivendo ad esempio nella 2 riga (relativa a 2 giorni) =$H$2/(1+($H$3*EXP(-$H$4*A2))). $H2$ significa che uso il parametro K inserito nella cella H2 etc.



A
B
C
D
E
F
G
H
1
giorni
popolazione
modello
err quadr
mod di cresc logistica
parametri:

2
2
25
23,8585963

y=K/(1+a exp(-bt))
K =
570
3
3
51
42,026982



a =
76
4
4
67
72,2015766



b =
0,6
5
5
98
119,15171





6
6
197
185,269278





7
7
287
266,396987





8
8
339
350,669797





9
9
409
424,340735





10
10
512
479,642473





11
11
519
516,590732





12
12
534
539,394476





13
13
544
552,786313





14
14
551
560,422427





15
15
556
564,703563






Successivamente calcoliamo l'errore quadratico inserendo in ogni riga della colonna D formule come  =(C2-B2)^2 (questo per la cella D2) etc.

Otteniamo:


A
B
C
D
E
F
G
H
1
giorni
popolazione
modello
err quadr
mod di cresc logistica
parametri:

2
2
25
23,8585963
1,30280247
y=K/(1+a exp(-bt))
K =
570
3
3
51
42,026982
80,5150525


a =
76
4
4
67
72,2015766
27,0563993


b =
0,6
5
5
98
119,15171
447,394843




6
6
197
185,269278
137,609832




7
7
287
266,396987
424,484139




8
8
339
350,669797
136,184158




9
9
409
424,340735
235,338149




10
10
512
479,642473
1047,00953




11
11
519
516,590732
5,8045703




12
12
534
539,394476
29,1003731




13
13
544
552,786313
77,1992874




14
14
551
560,422427
88,7821387




15
15
556
564,703563
77,1992874




16


sum err^2
2814,98057














Usando il "Solver" ("risolutore") di Excel (che si trova sotto il menu "Tools" - "Strumenti") ottimizziamo i parametri K, a e b minimizzando la somma degli scarti quadratici medi. Bisogna allora impostare il "Solver" indicacando che i parametri sono nelle celle e H2-H4 (e quindi inserendo il comando $H$2:$H$4 sotto "changing cells") degli scarti quadratici medie la somma degli scarti quadratici medi D16  ("Target cell" D16).

Si trova:



A
B
C
D
E
F
G
H
1
giorni
popolazione
modello
errore quadr
mod di cresc logistica
parametri:

2
2
25
21,6673197 11,10675795 y=K/(1+a exp(-bt))

K =
560,7788183
3
3
51
39,2849303
137,242858


a =
87,4125397
4
4
67
69,3837439
5,68223513


b =
0,62826073
5
5
98
117,353829
374,57069




6
6
197
185,939714
122,329929




7
7
287
270,185829
282,716347




8
8
339
356,31798
299,912445




9
9
409
429,340015
413,716198




10
10
512
482,04526
897,286451




11
11
519
515,829028
10,0550656




12
12
534
535,865722
3,48091893




13
13
544
547,205928
10,2779714




14
14
551
553,454731
6,02570553




15
15
556
556,847327
10,2779714




16


sum err^2
2584,68154






   
 

Esercizi:


  • Ripetere la procedura fatta per il modello sigmoide.
  • Analizzare in modo analogo 
(1) la crescita del Saccharomyces cerevisiae come nella seguente tabella (file excel)

tempo volume
0
1,5
9
10
18
18,5
23
25,5
27
34
38
42
45,5
47
0,37
1,63
6,2
8,87
10,66
10,97
12,5
12,6
12,9
13,27
12,77
12,87
12,9
12,7

(2) la crescita di girasoli come nella seguente tabella

tempo
cm
7.0
14.0
21.0
28.0
35.0
42.0
49.0
56.0
63.0
70.0
77.0
84.0
17.93
36.36
67.76
98.1
131.0
169.5
205.5
228.3
247.1
250.5
253.8
254.5

  • Determinare la retta di regressione per la crescita del bacillo aspergillus niger (file excel)
    Tempo 
    Diametro
    3
    4
    5
    6
    7
    8
    9
    9,3
    16,8
    22,8
    28,5
    33,6
    36,6
    42,8