Regressione non lineare in Python: 6 modelli con asintoto, linearizzazione e curve_fit

Non tutti i fenomeni crescono indefinitamente.

In molti problemi reali una variabile osservata cresce rapidamente, poi rallenta e tende progressivamente verso un valore limite.

Quel valore prende il nome di asintoto.

È un comportamento che incontriamo in moltissimi contesti:

• crescita di popolazioni;
• processi biologici;
• cinetica enzimatica;
• assorbimento e saturazione;
• crescita di microrganismi;
• diffusione di innovazioni;
• fenomeni industriali;
• processi di apprendimento;
• modelli farmacocinetici.

Dal punto di vista della data science, questi fenomeni richiedono spesso modelli di regressione non lineare, nei quali uno o più parametri compaiono in modo non lineare.

Tra i modelli più utilizzati troviamo:

• [math]y = A – B e^{-k x}[/math]   modello esponenziale saturante,
• [math]y = A – \frac{B}{x}[/math]   modello iperbolico,
• [math]v = \frac{V_{\max} S}{K_m + S}[/math]   modello di Michaelis-Menten,
• [math]y = \frac{K}{1 + e^{-(g x + d)}}[/math]   modello logistico,
• [math]y = A \left(1 – e^{-k x^n}\right)[/math]   modello di saturazione a potenza,
• e [math]V = A e^{-B e^{-c t}}[/math]   modello di Gompertz.

Alcuni di questi modelli possono essere trasformati in una relazione lineare. Ma questo non significa automaticamente che la regressione sulla variabile trasformata sia il modo migliore per stimare i parametri.

Questa distinzione è fondamentale.

La linearizzazione è uno strumento matematico. La regressione non lineare stima invece direttamente il modello sulla scala originale dei dati.

In questo articolo useremo le linearizzazioni soprattutto per capire la struttura matematica dei modelli e, quando opportuno, per ottenere stime iniziali. Il fit definitivo verrà effettuato, quando necessario, sul modello originale mediante scipy.optimize.curve_fit.

L’obiettivo non è quindi soltanto trovare dei numeri, ma capire perché un metodo di stima è preferibile a un altro.

1. Prima distinzione: modello non lineare e trasformazione lineare

Supponiamo di avere un modello

[math]y = A – B e^{-k x}.[/math]

Il modello è non lineare nei parametri perché [math]k[/math] compare nell’esponente.

Possiamo però scrivere:

[math]A – y = B e^{-k x}[/math]

e quindi:

[math]\ln(A – y) = \ln B – k x.[/math]

Sembra dunque che basti fare una regressione lineare.

Ma attenzione.

La regressione lineare sulla quantità trasformata minimizza:

[math]\displaystyle \sum_i \left[ \ln(A – y_i) – \ln(A – \hat y_i) \right]^2.[/math]

La regressione non lineare sul modello originale minimizza invece:

[math]\displaystyle \sum_i (y_i – \hat y_i)^2.[/math]

Le due funzioni obiettivo non sono equivalenti.

Questo significa che possiamo ottenere parametri diversi.

Regola generale

Se trasformiamo la variabile dipendente, cambiamo anche il modo in cui gli errori vengono pesati.

Per questo una linearizzazione può essere molto utile per:

• comprendere il modello;
• costruire un grafico diagnostico;
• ottenere valori iniziali;
• verificare qualitativamente la plausibilità del modello.

Ma non dobbiamo automaticamente considerarla il metodo finale di stima.

Esercizio 1 – Crescita batterica con asintoto

Testo

In un esperimento di laboratorio si misura il numero di colonie batteriche, espresso in migliaia, in funzione del tempo.

I dati sono:

Supponiamo che la crescita possa essere descritta da:

[math]y(t) = A – B e^{-k t}.[/math]

Determinare [math]A[/math], [math]B[/math] e [math]k[/math].

Soluzione

Passo 1 – Interpretare i parametri

Il modello è:

[math]y(t) = A – B e^{-k t}.[/math]

Se [math]t \to \infty[/math], allora [math]e^{-k t} \to 0[/math] e quindi [math]y(t) \to A.[/math] Pertanto [math]A[/math] rappresenta l’asintoto superiore.

Dai dati vediamo che le ultime osservazioni sono vicine a 10: [math]9.4,\;9.7,\;9.9,\;10.0.[/math] Possiamo quindi utilizzare [math]A_0 = 10[/math] come prima stima dell’asintoto.

Passo 2 – Linearizzazione

Partiamo da [math]y = A – B e^{-k t}[/math]. Portiamo [math]y[/math] dall’altra parte: [math]A – y = B e^{-k t}[/math]. Applicando il logaritmo: [math]\ln(A – y) = \ln B – k t.[/math] Questa è una retta nella variabile [math]t[/math]:

[math]z = \alpha + \beta t[/math] con [math]z = \ln(A – y)[/math], [math]\alpha = \ln B[/math], [math]\beta = -k[/math]. Quindi [math]B = e^{\alpha}[/math] e [math]k = -\beta.[/math]

Attenzione: per [math]t=10[/math] si ha [math]y=10[/math] e [math]A-y=0[/math]; il logaritmo non è definito. Per questa ragione utilizziamo la linearizzazione soltanto sui primi dieci punti.

Passo 3 – Calcolare la regressione linearizzata

Utilizzando [math]A_0=10[/math]:

La regressione lineare produce approssimativamente:

[math]\beta = -0.4624[/math] e [math]\alpha = 2.5073.[/math]

Quindi:

[math]k_0 = -\beta \approx 0.4624[/math] e [math]B_0 = e^{2.5073} \approx 12.27.[/math]

Otteniamo dunque le stime iniziali:

[math]A_0 = 10,\quad B_0 = 12.27,\quad k_0 = 0.4624.[/math]

Passo 4 – Perché non ci fermiamo qui?

Se usassimo direttamente queste stime avremmo:

[math]y(t) = 10 – 12.27 e^{-0.4624 t}.[/math]

Per [math]t=0[/math]: [math]y(0) = 10 – 12.27 = -2.27[/math]. Il modello prevederebbe quindi un numero negativo di colonie. Questo risultato è chiaramente incompatibile con i dati. Non significa che la matematica della trasformazione sia sbagliata. Significa che la regressione linearizzata sta ottimizzando una funzione obiettivo diversa da quella che vogliamo utilizzare sulla scala originale. La linearizzazione ci ha dato un buon punto di partenza concettuale, ma non necessariamente il miglior modello finale.

Passo 5 – Regressione non lineare sul modello originale

Utilizziamo curve_fit.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# -----------------------------
# Dati sperimentali
# -----------------------------
t = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float)
y = np.array([
    1.8, 3.3, 5.0, 6.4, 7.5,
    8.4, 9.0, 9.4, 9.7, 9.9, 10.0
])

# -----------------------------
# Modello esponenziale saturante
# -----------------------------
def modello(t, A, B, k):
    return A - B * np.exp(-k * t)

# -----------------------------
# Stime iniziali ottenute dalla linearizzazione
# -----------------------------
p0 = [10.0, 12.27, 0.4624]

# -----------------------------
# Regressione non lineare
# -----------------------------
popt, pcov = curve_fit(modello, t, y, p0=p0)

A, B, k = popt

print(f"A = {A:.4f}")
print(f"B = {B:.4f}")
print(f"k = {k:.4f}")
    

Il risultato è circa:

A = 11.0558
B = 9.5047
k = 0.2413
    

Il modello finale diventa quindi:

[math]\boxed{ y(t) \approx 11.056 – 9.505 e^{-0.241 t} }[/math]

Passo 6 – Controllare il modello

Per [math]t=0[/math]: [math]y(0) = 11.056 – 9.505 \approx 1.55[/math]. Il primo dato osservato è 1.8. Per [math]t \to \infty[/math]: [math]y(t) \to 11.056[/math]. L’asintoto stimato è quindi coerente con l’andamento dei dati, ma leggermente superiore al valore iniziale ipotizzato.

Passo 7 – Valutare quantitativamente il fit

Possiamo calcolare SSE e RMSE.

# Predizioni
y_pred = modello(t, A, B, k)

# Residui
residui = y - y_pred

# SSE
SSE = np.sum(residui**2)

# RMSE
RMSE = np.sqrt(np.mean(residui**2))

print(f"SSE = {SSE:.4f}")
print(f"RMSE = {RMSE:.4f}")
    

Le definizioni sono:

[math]\displaystyle SSE = \sum_i (y_i – \hat y_i)^2[/math]

e

[math]\displaystyle RMSE = \sqrt{\frac{SSE}{n}}.[/math]

L’SSE misura la quantità totale di errore quadratico. L’RMSE è più facilmente interpretabile perché mantiene la stessa unità di misura della variabile [math]y[/math].

Passo 8 – Analizzare i residui

plt.figure(figsize=(9, 5))
plt.axhline(0, linestyle="--")
plt.scatter(t, residui, s=70)
plt.xlabel("Tempo (ore)")
plt.ylabel("Residuo")
plt.title("Residui del modello esponenziale saturante")
plt.grid(True)
plt.show()
    

Il grafico dei residui è importante perché un modello non deve essere giudicato soltanto da quanto la curva “sembra passare vicino ai dati”. Se i residui mostrano una struttura sistematica, potrebbe esserci un problema di specificazione del modello.

Passo 9 – Visualizzare il fit

t_fit = np.linspace(0, 11, 300)
y_fit = modello(t_fit, A, B, k)

plt.figure(figsize=(10, 6))
plt.scatter(t, y, s=70, label="Dati sperimentali")
plt.plot(t_fit, y_fit, linewidth=3, label="Regressione non lineare")
plt.axhline(A, linestyle=":", linewidth=2, label=f"Asintoto A = {A:.2f}")
plt.xlabel("Tempo (ore)")
plt.ylabel("Colonie (migliaia)")
plt.title("Crescita batterica con asintoto")
plt.legend()
plt.grid(True)
plt.show()
    

Conclusione dell’esercizio

La linearizzazione è stata utile per ottenere [math]A_0, B_0, k_0[/math]. Ma la stima finale è stata effettuata direttamente sul modello originale.

Questo è il principio che utilizzeremo anche negli esercizi successivi.

Esercizio 2 – Modello iperbolico

Testo

In un processo di cristallizzazione la velocità di accrescimento [math]v[/math], espressa in µm/min, dipende dal tempo [math]t[/math] secondo:

[math]v(t) = A – \frac{B}{t}.[/math]

I dati sono:

Stimare [math]A[/math] e [math]B[/math].

Soluzione

Passo 1 – Una linearizzazione particolare

Definiamo [math]u = \frac{1}{t}[/math]. Il modello diventa [math]v = A – B u[/math]. Possiamo quindi scrivere [math]v = A + C u[/math] dove [math]C = -B[/math]. Questa volta la situazione è diversa dall’Esercizio 1. La trasformazione riguarda la variabile indipendente, non la variabile dipendente. La regressione rimane quindi una regressione lineare di [math]v[/math] su [math]u[/math].

Passo 2 – Applicare OLS

import numpy as np
import matplotlib.pyplot as plt

t = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float)
v = np.array([2.0, 4.5, 5.7, 6.3, 6.7, 7.0, 7.2, 7.3, 7.4, 7.5])

# Trasformazione della variabile indipendente
u = 1 / t

# Regressione lineare: v = A + C*u
C, A = np.polyfit(u, v, 1)

# Nella formulazione originale:
B = -C

print(f"A = {A:.4f}")
print(f"B = {B:.4f}")
print(f"C = {C:.4f}")
    

Otteniamo:

[math]A \approx 7.9793[/math] e [math]B \approx 6.2113.[/math]

Il modello è quindi:

[math]\boxed{ v(t) \approx 7.979 – \frac{6.211}{t} }[/math]

Passo 3 – Interpretare l’asintoto

Quando [math]t \to \infty[/math], [math]\frac{B}{t} \to 0[/math], quindi [math]v(t) \to A[/math]. L’asintoto è:

[math]\boxed{A \approx 7.98\ \mu m/min}[/math]

Passo 4 – Verificare alcune predizioni

Per [math]t=1[/math]: [math]v(1) = 7.979 – 6.211 \approx 1.768[/math] (dato osservato 2.0). Per [math]t=10[/math]: [math]v(10) = 7.979 – \frac{6.211}{10} \approx 7.358[/math] (dato osservato 7.5). Il modello riproduce quindi ragionevolmente l’andamento.

Passo 5 – Grafico

t_fit = np.linspace(0.5, 11, 300)
v_fit = A - B / t_fit

plt.figure(figsize=(10, 6))
plt.scatter(t, v, s=70, label="Dati sperimentali")
plt.plot(t_fit, v_fit, linewidth=3, label=f"v = {A:.2f} - {B:.2f}/t")
plt.axhline(A, linestyle="--", linewidth=2, label=f"Asintoto A = {A:.2f}")
plt.xlabel("Tempo (min)")
plt.ylabel("Velocità (µm/min)")
plt.title("Modello iperbolico con asintoto")
plt.xlim(0, 11)
plt.ylim(0, 9)
plt.legend()
plt.grid(True)
plt.show()
    

Una differenza importante rispetto all’Esercizio 1: Qui la trasformazione [math]u = 1/t[/math] non trasforma la variabile risposta [math]v[/math]. Per questo non introduciamo lo stesso tipo di problema incontrato con il logaritmo dell’Esercizio 1.

Non tutte le linearizzazioni hanno le stesse conseguenze statistiche.

Esercizio 3 – Michaelis-Menten

Testo

La velocità di una reazione enzimatica dipende dalla concentrazione del substrato secondo:

[math]v = \frac{V_{\max} S}{K_m + S}.[/math]

I dati sono:

Stimare [math]V_{\max}[/math] e [math]K_m[/math].

Soluzione

Passo 1 – Lineweaver-Burk

Partiamo da [math]v = \frac{V_{\max} S}{K_m + S}[/math]. Prendiamo il reciproco:

[math]\displaystyle \frac{1}{v} = \frac{K_m + S}{V_{\max} S} = \frac{K_m}{V_{\max}} \frac{1}{S} + \frac{1}{V_{\max}}.[/math]

Poniamo [math]x = \frac{1}{S}[/math] e [math]y = \frac{1}{v}[/math]. Otteniamo una retta [math]y = a + b x[/math] con [math]a = \frac{1}{V_{\max}}[/math] e [math]b = \frac{K_m}{V_{\max}}.[/math]

Passo 2 – Calcolare le stime iniziali

import numpy as np

S = np.array([0.1, 0.2, 0.5, 1, 2, 5, 10], dtype=float)
v = np.array([0.5, 0.8, 1.2, 1.6, 2.0, 2.4, 2.6])

x = 1 / S
y = 1 / v

pendenza, intercetta = np.polyfit(x, y, 1)   # y = intercetta + pendenza*x

Vmax_0 = 1 / intercetta
Km_0 = pendenza * Vmax_0

print(f"Vmax iniziale = {Vmax_0:.4f}")
print(f"Km iniziale = {Km_0:.4f}")
    

Le stime iniziali sono circa:

[math]V_{\max,0} \approx 2.33[/math] e [math]K_{m,0} \approx 0.37.[/math]

Passo 3 – Perché non fermarsi alla linearizzazione?

Il problema della trasformazione di Lineweaver-Burk è che stiamo trasformando la variabile dipendente: [math]v \to \frac{1}{v}[/math]. Gli errori non vengono più trattati nello stesso modo. In particolare, le osservazioni con [math]v[/math] piccolo possono acquisire un peso molto maggiore dopo il reciproco. Quindi la retta può essere utile per ottenere valori iniziali, ma per la stima finale è preferibile lavorare sul modello originale.

Passo 4 – Regressione non lineare

from scipy.optimize import curve_fit

def michaelis_menten(S, Vmax, Km):
    return Vmax * S / (Km + S)

p0 = [Vmax_0, Km_0]
popt, pcov = curve_fit(
    michaelis_menten,
    S,
    v,
    p0=p0,
    bounds=(0, np.inf)
)

Vmax, Km = popt
print(f"Vmax = {Vmax:.4f}")
print(f"Km = {Km:.4f}")
    

Qui imponiamo [math]V_{\max} > 0[/math] e [math]K_m > 0[/math]. L’uso di bounds è una buona pratica perché incorpora nel problema numerico una conoscenza fisica del modello.

Passo 5 – Valutare il modello

v_pred = michaelis_menten(S, Vmax, Km)
residui = v - v_pred
SSE = np.sum(residui**2)
RMSE = np.sqrt(np.mean(residui**2))
print(f"SSE = {SSE:.4f}")
print(f"RMSE = {RMSE:.4f}")
    

Passo 6 – Grafico

S_fit = np.linspace(0.05, 11, 300)
v_fit = michaelis_menten(S_fit, Vmax, Km)

plt.figure(figsize=(10, 6))
plt.scatter(S, v, s=70, label="Dati sperimentali")
plt.plot(S_fit, v_fit, linewidth=3, label="Michaelis-Menten")
plt.axhline(Vmax, linestyle="--", linewidth=2, label=f"Vmax = {Vmax:.2f}")
plt.xlabel("Concentrazione del substrato S (mM)")
plt.ylabel("Velocità v")
plt.title("Fit non lineare del modello di Michaelis-Menten")
plt.legend()
plt.grid(True)
plt.show()
    

Interpretazione

[math]V_{\max}[/math] rappresenta la velocità limite raggiunta quando il substrato diventa molto abbondante. [math]K_m[/math] determina la scala di concentrazione alla quale il sistema raggiunge una parte significativa della velocità massima.

Il punto fondamentale dal punto di vista statistico è però un altro:

La trasformazione Lineweaver-Burk può essere utilizzata per comprendere il modello e ottenere valori iniziali, mentre la stima finale viene effettuata sulla scala originale.

Esercizio 4 – Modello logistico

Testo

Una popolazione batterica, espressa in migliaia, segue approssimativamente il modello:

[math]y(t) = \frac{K}{1 + e^{-(g t + d)}}.[/math]

I dati sono:

Stimare [math]K, g, d[/math].

Soluzione

Passo 1 – Il problema della trasformazione logit

Dal modello [math]y = \frac{K}{1 + e^{-(g t + d)}}[/math] otteniamo [math]\frac{K – y}{y} = e^{-(g t + d)}[/math] e quindi [math]\ln\left( \frac{K – y}{y} \right) = -d – g t.[/math] Questa trasformazione è lineare solo se [math]K[/math] è già noto. Nel nostro caso [math]K[/math] è uno dei parametri da stimare. Non possiamo quindi affermare rigorosamente di aver stimato contemporaneamente [math]K, g, d[/math] con una semplice regressione lineare.

Passo 2 – Utilizzare [math]K=10[/math] come valore iniziale

Dai dati: [math]y \to 10[/math]. Possiamo quindi utilizzare [math]K_0 = 10[/math] come stima iniziale. Costruiamo [math]z = \ln\left( \frac{K_0 – y}{y} \right).[/math] Attenzione: per [math]t=9[/math] e [math]t=10[/math] si ha [math]y=10.0[/math]. Con [math]K_0=10[/math] il rapporto diventa zero e il logaritmo non è definito. Per questo escludiamo gli ultimi due punti dalla regressione linearizzata, esattamente come abbiamo fatto nell’Esercizio 1.

t = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float)
y = np.array([0.5, 1.2, 2.5, 4.1, 6.0, 7.8, 9.0, 9.6, 9.9, 10.0, 10.0])

K0 = 10.0

# Escludiamo gli ultimi due punti per evitare log(0)
t_lin = t[:-2]
y_lin = y[:-2]
z = np.log((K0 - y_lin) / y_lin)

b, a = np.polyfit(t_lin, z, 1)   # z = a + b*t

d0 = -a
g0 = -b

print(f"K iniziale = {K0:.2f}")
print(f"g iniziale = {g0:.4f}")
print(f"d iniziale = {d0:.4f}")
    

Queste stime servono principalmente a inizializzare l’ottimizzazione.

Passo 3 – Fit non lineare

Definiamo direttamente il modello:

def logistico(t, K, g, d):
    return K / (1 + np.exp(-(g * t + d)))
    

Ora possiamo stimare contemporaneamente i tre parametri.

p0 = [10.0, g0, d0]
popt, pcov = curve_fit(
    logistico,
    t,
    y,
    p0=p0,
    bounds=([0, 0, -np.inf], [np.inf, np.inf, np.inf])
)

K, g, d = popt
print(f"K = {K:.4f}")
print(f"g = {g:.4f}")
print(f"d = {d:.4f}")
    

La differenza concettuale è fondamentale:

• Metodo linearizzato: [math]K[/math] viene fissato.
• Regressione non lineare: [math]K, g, d[/math] vengono stimati contemporaneamente.

Passo 4 – Valutare il modello

y_pred = logistico(t, K, g, d)
residui = y - y_pred
SSE = np.sum(residui**2)
RMSE = np.sqrt(np.mean(residui**2))
print(f"SSE = {SSE:.4f}")
print(f"RMSE = {RMSE:.4f}")
    

Passo 5 – Grafico

t_fit = np.linspace(0, 10, 300)
y_fit = logistico(t_fit, K, g, d)

plt.figure(figsize=(10, 6))
plt.scatter(t, y, s=70, label="Dati")
plt.plot(t_fit, y_fit, linewidth=3, label="Fit logistico")
plt.axhline(K, linestyle="--", linewidth=2, label=f"Asintoto K = {K:.2f}")
plt.xlabel("Tempo")
plt.ylabel("Popolazione (migliaia)")
plt.title("Modello logistico")
plt.legend()
plt.grid(True)
plt.show()
    

Interpretazione

Il modello logistico descrive tre fasi: crescita iniziale relativamente lenta, crescita rapida, rallentamento verso [math]K[/math]. Il parametro [math]K[/math] rappresenta la capacità portante o il livello massimo del modello.

Esercizio 5 – Modello di saturazione a potenza

Testo

Consideriamo il modello:

[math]C(t) = A \left(1 – e^{-k t^n}\right).[/math]

I dati sono:

Stimare [math]A, k, n[/math].

Soluzione

Passo 1 – Linearizzazione

Partiamo da [math]C = A \left(1 – e^{-k t^n}\right)[/math]. Dividiamo per [math]A[/math]: [math]\frac{C}{A} = 1 – e^{-k t^n}[/math], quindi [math]1 – \frac{C}{A} = e^{-k t^n}[/math]. Prendendo il logaritmo: [math]\ln\left(1 – \frac{C}{A}\right) = -k t^n[/math]; cambiando segno: [math]-\ln\left(1 – \frac{C}{A}\right) = k t^n[/math]; prendendo nuovamente il logaritmo: [math]\ln\left[-\ln\left(1 – \frac{C}{A}\right)\right] = \ln k + n \ln t.[/math] Abbiamo quindi una relazione lineare [math]y = a + n x[/math] con [math]x = \ln t[/math] e [math]y = \ln\left[-\ln\left(1 – \frac{C}{A}\right)\right].[/math]

Passo 2 – Il problema dell’asintoto

Ancora una volta [math]A[/math] è sconosciuto. Possiamo osservare che i dati arrivano a 6.4, quindi un valore iniziale ragionevole può essere [math]A_0 = 6.8[/math]. Ma non dobbiamo interpretare 6.8 come una stima definitiva. È soltanto un valore iniziale.

Passo 3 – Calcolare le stime iniziali

t = np.array([0.5, 1, 1.5, 2, 2.5, 3, 4], dtype=float)
C = np.array([0.8, 2.5, 4.0, 5.0, 5.6, 6.0, 6.4])

A0 = 6.8
x = np.log(t)
y_trans = np.log(-np.log(1 - C / A0))

n0, log_k0 = np.polyfit(x, y_trans, 1)
k0 = np.exp(log_k0)

print(f"A iniziale = {A0:.2f}")
print(f"k iniziale = {k0:.4f}")
print(f"n iniziale = {n0:.4f}")
    

Questi valori non rappresentano ancora la soluzione finale.

Passo 4 – Regressione non lineare

Definiamo il modello originale:

def potenza_saturazione(t, A, k, n):
    return A * (1 - np.exp(-k * t**n))
    

Ora stimiamo direttamente [math]A, k, n[/math].

p0 = [A0, k0, n0]
popt, pcov = curve_fit(
    potenza_saturazione,
    t,
    C,
    p0=p0,
    bounds=(0, np.inf)
)

A, k, n = popt
print(f"A = {A:.4f}")
print(f"k = {k:.4f}")
print(f"n = {n:.4f}")
    

In questo modo l’asintoto non viene più imposto arbitrariamente. Viene stimato dal modello insieme agli altri parametri.

Passo 5 – Metriche

C_pred = potenza_saturazione(t, A, k, n)
residui = C - C_pred
SSE = np.sum(residui**2)
RMSE = np.sqrt(np.mean(residui**2))
print(f"SSE = {SSE:.4f}")
print(f"RMSE = {RMSE:.4f}")
    

Passo 6 – Grafico

t_fit = np.linspace(0.1, 4.5, 300)
C_fit = potenza_saturazione(t_fit, A, k, n)

plt.figure(figsize=(10, 6))
plt.scatter(t, C, s=70, label="Dati sperimentali")
plt.plot(t_fit, C_fit, linewidth=3, label="Fit non lineare")
plt.axhline(A, linestyle="--", linewidth=2, label=f"Asintoto A = {A:.2f}")
plt.xlabel("Tempo")
plt.ylabel("Concentrazione")
plt.title("Modello di saturazione a potenza")
plt.legend()
plt.grid(True)
plt.show()
    

Cosa abbiamo imparato?

Questo esercizio mostra chiaramente perché la stima di [math]A[/math] è delicata. Se [math]A[/math] viene scelto troppo basso, può accadere che [math]C > A[/math]. In quel caso [math]1 – C/A < 0[/math] e la trasformazione logaritmica non è più definita. La regressione non lineare evita questo problema concettuale perché lavora direttamente sul modello originale.

Esercizio 6 – Modello di Gompertz

Nota importante sui dati

L’esempio utilizza la crescita di un tumore soltanto come scenario matematico.

I dati sono costruiti a fini didattici e non rappresentano osservazioni cliniche reali. Il modello non deve essere interpretato come strumento diagnostico o prognostico.

Testo

Consideriamo:

[math]V(t) = A e^{-B e^{-c t}}.[/math]

I dati didattici sono:

Stimare [math]A, B, c[/math].

Soluzione

Passo 1 – Interpretare il modello

Il modello è [math]V(t) = A e^{-B e^{-c t}}[/math]. Quando [math]t \to \infty[/math], [math]e^{-c t} \to 0[/math], quindi [math]V(t) \to A e^0 = A[/math]. Pertanto [math]A[/math] è l’asintoto superiore.

Passo 2 – Linearizzazione

Dividiamo per [math]A[/math]: [math]\frac{V}{A} = e^{-B e^{-c t}}[/math]. Prendiamo il logaritmo: [math]\ln\left( \frac{V}{A} \right) = -B e^{-c t}[/math]; cambiamo segno: [math]-\ln\left( \frac{V}{A} \right) = B e^{-c t}[/math]; prendiamo ancora il logaritmo: [math]\ln\left[ -\ln\left( \frac{V}{A} \right) \right] = \ln B – c t.[/math] Definiamo [math]z = \ln\left[ -\ln\left( \frac{V}{A} \right) \right][/math]. Otteniamo [math]z = \ln B – c t[/math], che è una retta.

Passo 3 – Stima iniziale di [math]A[/math]

Dai dati osserviamo [math]V(30) = 9.3[/math]. Un valore iniziale ragionevole è [math]A_0 = 10[/math]. Ancora una volta: 10 non è la stima definitiva di [math]A[/math]: è soltanto un valore iniziale.

Passo 4 – Calcolare la trasformazione

t = np.array([0, 5, 10, 15, 20, 25, 30], dtype=float)
V = np.array([0.5, 1.2, 2.8, 5.0, 7.0, 8.5, 9.3])

A0 = 10.0
z = np.log(-np.log(V / A0))

b, a = np.polyfit(t, z, 1)   # z = a + b*t

c0 = -b
B0 = np.exp(a)

print(f"A iniziale = {A0:.2f}")
print(f"B iniziale = {B0:.4f}")
print(f"c iniziale = {c0:.4f}")
    

Queste stime permettono di costruire un punto iniziale per l’ottimizzazione.

Passo 5 – Fit non lineare

Definiamo il modello originale:

def gompertz(t, A, B, c):
    return A * np.exp(-B * np.exp(-c * t))
    

Ora stimiamo simultaneamente i tre parametri:

p0 = [A0, B0, c0]
popt, pcov = curve_fit(
    gompertz,
    t,
    V,
    p0=p0,
    bounds=(0, np.inf)
)

A, B, c = popt
print(f"A = {A:.4f}")
print(f"B = {B:.4f}")
print(f"c = {c:.4f}")
    

La differenza rispetto alla linearizzazione è importante: [math]A[/math] non viene più fissato a 10. L’algoritmo cerca simultaneamente [math]A, B, c[/math] che minimizzano l’errore quadratico sui dati originali.

Passo 6 – Calcolare SSE e RMSE

V_pred = gompertz(t, A, B, c)
residui = V - V_pred
SSE = np.sum(residui**2)
RMSE = np.sqrt(np.mean(residui**2))
print(f"SSE = {SSE:.4f}")
print(f"RMSE = {RMSE:.4f}")
    

Passo 7 – Grafico

t_fit = np.linspace(0, 35, 300)
V_fit = gompertz(t_fit, A, B, c)

plt.figure(figsize=(10, 6))
plt.scatter(t, V, s=70, label="Dati didattici")
plt.plot(t_fit, V_fit, linewidth=3, label="Fit Gompertz")
plt.axhline(A, linestyle="--", linewidth=2, label=f"Asintoto A = {A:.2f}")
plt.xlabel("Tempo (giorni)")
plt.ylabel("Volume (cm³)")
plt.title("Modello di Gompertz")
plt.legend()
plt.grid(True)
plt.show()
    

Come valutare una regressione non lineare

Trovare i parametri non è sufficiente.

Un modello può produrre una curva apparentemente convincente e tuttavia essere statisticamente inadeguato.

Per questo dobbiamo guardare almeno tre aspetti: errore complessivo, residui e plausibilità dei parametri.

SSE

La somma dei quadrati degli errori è:

[math]\displaystyle SSE = \sum_{i=1}^n (y_i – \hat y_i)^2.[/math]

Più è piccola, meglio il modello riproduce i dati, a parità di struttura del problema. Non è però una metrica sufficiente per confrontare modelli con caratteristiche molto diverse.

RMSE

L’errore quadratico medio può essere calcolato in due modi.

La versione più comune in data science è:

[math]\displaystyle RMSE = \sqrt{ \frac{1}{n} \sum_{i=1}^n (y_i – \hat y_i)^2 }.[/math]

In statistica inferenziale, però, si preferisce correggere per i gradi di libertà:

[math]\displaystyle RMSE = \sqrt{\frac{SSE}{n-p}}[/math]

dove [math]p[/math] è il numero di parametri stimati. Questa seconda versione è più appropriata quando si confrontano modelli con diverso numero di parametri.

Il vantaggio dell’RMSE è che è espresso nella stessa unità di misura di [math]y[/math]. Se, per esempio, [math]RMSE = 0.20[/math] e [math]y[/math] è espresso in cm³, possiamo interpretarlo come un errore tipico dell’ordine di 0.20 cm³, con le dovute cautele statistiche.

Residui

Definiamo [math]e_i = y_i – \hat y_i[/math]. Il grafico dei residui è fondamentale.

residui = y - y_pred
plt.figure(figsize=(9, 5))
plt.scatter(y_pred, residui, s=70)
plt.axhline(0, linestyle="--")
plt.xlabel("Valori predetti")
plt.ylabel("Residui")
plt.title("Residui vs valori predetti")
plt.grid(True)
plt.show()
    

Idealmente non dovremmo osservare: una curva, una struttura crescente, una struttura decrescente, una forma a imbuto o gruppi sistematici.

Se compare una struttura, il problema potrebbe non essere il metodo di ottimizzazione. Potrebbe essere il modello matematico scelto.

popt, pcov = curve_fit(modello, x, y, p0=p0)

se = np.sqrt(np.diag(pcov))

for nome, valore, errore in zip(["A", "B", "k"], popt, se):
    print(
        f"{nome} = {valore:.4f} "
        f"+/- {1.96 * errore:.4f}"
    )
    

std = np.sqrt(np.diag(pcov))

corr = pcov / np.outer(std, std)

print("Matrice di correlazione:")
print(corr)
    

Matrice di correlazione:
[[ 1.         -0.62024802 -0.91870327]
 [-0.62024802  1.          0.85115826]
 [-0.91870327  0.85115826  1.        ]]

popt, pcov = curve_fit(...)
se = np.sqrt(np.diag(pcov))

Un modello può adattarsi bene ma essere comunque problematico

Supponiamo di ottenere un RMSE molto basso. Possiamo concludere che il modello è buono?

Non necessariamente.

Dobbiamo chiederci:

• I parametri hanno senso? Se il modello richiede [math]A > 0, k > 0[/math], possiamo imporre questi vincoli.
• I residui sono casuali? Se mostrano una struttura, il modello potrebbe essere mal specificato.
• I parametri sono identificabili? Con pochi dati è possibile che diversi valori dei parametri producano curve molto simili.
• Il modello generalizza? Un fit eccellente sui dati osservati non garantisce necessariamente buone previsioni fuori dall’intervallo osservato.

Linearizzazione o regressione non lineare?

Possiamo riassumere le situazioni incontrate negli esercizi.

La regola più importante è:

Non scegliere una trasformazione soltanto perché rende l’equazione più facile da calcolare. Chiediti sempre quale funzione di errore stai minimizzando dopo la trasformazione.

Uno schema Python riutilizzabile

Molti problemi di regressione non lineare possono essere affrontati con uno schema simile.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 1. Dati
x = np.array([...], dtype=float)
y = np.array([...], dtype=float)

# 2. Definizione del modello
def modello(x, p1, p2, p3):
    return ...

# 3. Stime iniziali
p0 = [p1_0, p2_0, p3_0]

# 4. Fit non lineare
popt, pcov = curve_fit(modello, x, y, p0=p0)

# 5. Parametri stimati
print(popt)

# 6. Predizioni
y_pred = modello(x, *popt)

# 7. Residui
residui = y - y_pred

# 8. Metriche
SSE = np.sum(residui**2)
RMSE = np.sqrt(np.mean(residui**2))
print("SSE:", SSE)
print("RMSE:", RMSE)

# 9. Errore standard dei parametri
se = np.sqrt(np.diag(pcov))
print("SE:", se)

# 10. Grafico
x_fit = np.linspace(x.min(), x.max(), 300)
y_fit = modello(x_fit, *popt)

plt.scatter(x, y, label="Dati")
plt.plot(x_fit, y_fit, label="Modello")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid(True)
plt.show()
    

Questo schema separa chiaramente:

[math]\boxed{ \text{modello} \to \text{stima} \to \text{predizione} \to \text{diagnostica} }[/math]

Cosa abbiamo imparato dai sei modelli

I sei esercizi mostrano sei situazioni diverse.

• Esponenziale saturante: [math]y = A – B e^{-k x}[/math] – La linearizzazione è utile per ottenere valori iniziali, ma la stima finale deve essere effettuata direttamente sul modello originale. I parametri ottenuti dalla linearizzazione possono portare a valori fisicamente implausibili (come [math]y(0) < 0[/math]).
• Iperbolico: [math]v = A – \frac{B}{t}[/math] – La trasformazione [math]u = 1/t[/math] produce una regressione lineare nei parametri senza trasformare la variabile risposta.
• Michaelis-Menten: [math]v = \frac{V_{\max} S}{K_m + S}[/math] – La trasformazione di Lineweaver-Burk è elegante matematicamente, ma può alterare fortemente il comportamento degli errori. Le stime iniziali ottenute vanno verificate.
• Logistico: [math]y = \frac{K}{1 + e^{-(g t + d)}}[/math] – La trasformazione logit richiede la conoscenza di [math]K[/math]. Se [math]K[/math] deve essere stimato, la regressione non lineare è una soluzione più naturale. Va prestata attenzione ai punti in cui [math]y[/math] raggiunge il valore iniziale di [math]K[/math].
• Saturazione a potenza: [math]C = A \left(1 – e^{-k t^n}\right)[/math] – La doppia trasformazione logaritmica dipende dall’asintoto [math]A[/math], che non è noto a priori.
• Gompertz: [math]V = A e^{-B e^{-c t}}[/math] – Anche qui la linearizzazione dipende dalla conoscenza di [math]A[/math], rendendo particolarmente interessante il confronto con il fit non lineare.

La lezione più importante

La regressione non lineare non consiste semplicemente nel trovare una curva che passa vicino ai punti. È un problema di stima statistica dei parametri di un modello matematico.

La linearizzazione può essere estremamente utile perché permette di:

• comprendere la struttura del modello;
• costruire diagnostiche;
• ottenere stime iniziali;
• semplificare alcuni problemi.

Ma bisogna sempre ricordare che una trasformazione può modificare la struttura degli errori e quindi la funzione obiettivo.

Per questo, quando possibile, è spesso preferibile stimare direttamente il modello originale:

[math]\displaystyle \boxed{ \min_{\theta} \sum_i \left[ y_i – f(x_i; \theta) \right]^2 }[/math]

dove [math]\theta[/math] rappresenta il vettore dei parametri.

In Python, uno strumento semplice e potente per questo scopo è:

scipy.optimize.curve_fit

Conclusione

Le curve con asintoto sono un ottimo esempio di come matematica, statistica e programmazione debbano lavorare insieme.

La matematica ci permette di costruire il modello. La statistica ci permette di stimarne i parametri e di quantificare l’incertezza. Python ci permette di automatizzare l’ottimizzazione, visualizzare i risultati e diagnosticare il comportamento del modello.

La domanda corretta, quindi, non è soltanto:

“Come posso trasformare questa curva in una retta?”

La domanda più importante è:

“Quale modello descrive il fenomeno e qual è il modo più appropriato per stimarne i parametri?”

Ed è proprio questa la differenza tra una semplice linearizzazione matematica e una vera analisi di data science.