Monte Carlo: dare una forma all'incertezza

Un sistema di ODE con parametri fissati restituisce una sola traiettoria. È rassicurante, ed è quasi sempre falso: ogni parametro arriva dalla letteratura con il suo intervallo, e spesso con un intervallo largo.

L’idea

Invece di un valore per parametro, si usa una distribuzione. Poi:

1. estraggo un set di parametri dalle distribuzioni;
2. integro il modello;
3. ripeto N volte;
4. guardo la distribuzione delle traiettorie, non una traiettoria.

Il risultato è una banda: la mediana e, per esempio, l’intervallo fra il 5° e il 95° percentile. Se due scenari hanno bande che si sovrappongono quasi del tutto, il modello non sa distinguerli — e devo dirlo.

# schema generale, non il codice di un progetto specifico
rng = np.random.default_rng(seed)
runs = []
for _ in range(N):
    theta = sample_parameters(rng)      # una estrazione dalle distribuzioni
    runs.append(simulate(theta, protocol))
runs = np.array(runs)
p5, p50, p95 = np.percentile(runs, [5, 50, 95], axis=0)

Scelte che contano più del codice

• La forma della distribuzione. Per parametri positivi con incertezza moltiplicativa (costanti di velocità, emivite) una log-normale è di solito più sensata di una normale, che ammetterebbe valori negativi.
• Le correlazioni. Due parametri stimati dallo stesso esperimento non sono indipendenti; trattarli come tali gonfia o sgonfia la banda.
• Il seed. Va fissato, altrimenti i risultati non sono riproducibili.
• N. Abbastanza grande da rendere stabili i percentili estremi. Nel modello mitocondriale la matrice di scenari usa N=400.

Cosa non fa

Monte Carlo propaga l’incertezza dei parametri, non quella della struttura. Se il modello ha sbagliato la forma di un’equazione, mille estrazioni non lo correggeranno. Per questo la banda va letta come un limite inferiore dell’incertezza vera — una distinzione di So, stimo, speculo, non so.

Lo stesso approccio sta dietro alle bande di Caffeine Lab, dove la variabilità individuale nel metabolismo della caffeina è enorme.