2  Recap — analisi dati

Qui puntiamo a ricapitolare alcune basi partendo dalle note Dispense di Analisi dati di Lab1, aggiungendo solo qualche dettaglio sulla matrice di covarianza e qualche esempio pertinente in più.

2.1 Precisione e accuratezza

Lo scopo primario di ogni misura è quantificare un parametro fisico, ossia produrre una stima accompagnata da una quantificazione del suo grado di incertezza. Secondo il Vocabolario Internazionale di Metrologia (JCGM 2012) ogni misura è caratterizzata da:

1. Precisione \(\implies\) grado di concordanza tra misurazioni ripetute secondo una data procedura. Incertezza di natura statistica: può essere quantificata dai dati e abbattuta acquisendo più dati.
2. Accuratezza \(\implies\) grado di concordanza tra valore misurato e valore vero del misurando. Incertezza non statistica, bensì sistematica: può essere quantificata confrontando diversi strumenti/metodi di misura, non cambia acquisendo più dati nello stesso modo.

import { renderAccuracy } from "./_interactive/ch02-accuracy.js"
{
  const container = document.createElement("div");
  container.style.cssText = "width:100%;position:relative;margin:0 auto;";
  yield container;
  renderAccuracy(container);
}

Figura 2.1: Classica rappresentazione grafica di precisione versus accuratezza.

DefinizioneCome riportare le incertezze

Restano valide e cruciali le regole già note, a cui è obbligatorio attenersi:

• Le incertezze si riportano con una o al massimo due cifre significative: andare oltre non ha senso;
• Il valore misurato si riporta con un numero di cifre significative consistente con l’incertezza.
• Vanno indicate sempre le unità.

Esempio: \(V_{out} = 10.6 \pm 0.4\,\mathrm{V}\). Le incertezze statistiche e sistematiche possono (non richiesto) essere riportate separatamente, nel qual caso si usa un doppio \(\pm A \pm B\) indicando prima l’errore statistico e poi quello sistematico.

AttenzioneMisurare di più: sì, no, quanto?

Sì e no. In questo laboratorio sarà facile fare molte (a volte anche troppe) misure e medie e la risposta è che spesso non è bene/utile misurare troppo. Come noto, l’incertezza statistica viene tendenzialmente abbattuta facendo delle medie, ma questo non abbatte gli errori sistematici e quindi ad un certo punto non migliora le vostre stime. Inoltre, una misura dominata da errori sistematici non è quasi mai desiderabile, dato che gli errori sistematici sono spesso di difficile gestione. Quindi il consiglio è di essere moderati: è vero che facendo molte misure una stima può migliorare, ma se bastasse questo per migliorare ad libitum non esisterebbero gli istituti di metrologia. Più realisticamente ad un certo punto finirete per trovarvi in situazioni di difficile gestione, l’argomento verrà ripreso nel contesto dei best fit.

Quanto migliora? Il punto di partenza è un risultato già noto: la media \(\mu\) di enne misure \(x_1, \cdots, x_n\) ripetute fluttua meno della singola misura, e precisamente

\[\sigma_\mu^2 = \mathrm{Var}\!\left(\frac{1}{n}\sum_{i=1}^n x_i\right) = \frac{1}{n^2}\sum_i \sigma^2 = \frac{\sigma^2}{n} \qquad\Longrightarrow\qquad \sigma_\mu = \frac{\sigma}{\sqrt n}.\]

Senza voler infrangere troppe certezze, è importante ricordare che questo vale solo se la varianza della somma è uguale alla somma delle varianze, il che richiede che gli errori siano statisticamente indipendenti. Dato che le misure sono tipicamente ripetute nel tempo, questo si traduce nella richiesta che il rumore non sia correlato nel tempo: vedremo che questo equivale a richiedere che il suo spettro non abbia una dipendenza dalla frequenza. Questa condizione non è affatto scontata e se fallisce mediare può ridurre le fluttuazioni più o meno velocemente di \(1/\sqrt n\) a seconda dei casi. Come esempio classico, Figura 2.2 riporta un segnale con rumore casuale più un disturbo a \(50\,\mathrm{Hz}\) con 100 campioni per periodo. Al variare del numero di acquisizioni la deviazione standard della media ha dei minimi o massimi per acquisizioni che durano un numero intero o semi-intero di periodi.

Per chiarire quanto sia semplice incappare in una correlazione temporale, si consideri che un semplice filtro passa-basso \(RC\) può facilmente generare segnali correlati nel tempo. Immaginiamoci che all’ingresso del filtro si manifesti una fluttuazione di voltaggio, che per essere scorrelata nel tempo deve essere sostanzialmente una \(\delta(t)\). In un filtro RC questo impulso genera un disturbo in uscita che decade nel tempo come \(e^{-t/\tau}\) con \(\tau=RC\): questa è una correlazione temporale, ed è già sufficiente ad intaccare la validità della legge \(1/\sqrt{n}\).

Per rendere le cose ancora meno ovvie, si consideri anche che nel caso in cui siano presenti dei drift o un rumore \(1/f\), mediare può dare un miglioramento inferiore a quello previsto dalla legge \(1/\sqrt{n}\) o perfino nessun miglioramento. Come messaggio finale: attendersi un abbattimento delle fluttuazioni statistiche secondo la legge \(1/\sqrt{n}\) non è irragionevole, ma è sempre meglio verificarlo e non crederci senza aver compreso davvero il segnale che si sta misurando.

import { renderMediare } from "./_interactive/ch02-mediare.js"
{
  const container = document.createElement("div");
  container.style.cssText = "width:100%;position:relative;margin:0 auto;";
  yield container;
  renderMediare(container);
}

Figura 2.2: Segnale costante in presenza di rumore correlato nel tempo. Trascinando la zona evidenziata si cambia il numero di punti su cui si fa la media e a destra si vede l’istogramma delle misure ripetute.

Riassumiamo rapidamente a seguire altri strumenti e concetti utili.

Definizione 2.1: Distribuzioni, media e varianza

Per una variabile aleatoria \(x\) con densità di probabilità \(P(x)\), il valore di aspettazione di una funzione \(f(x)\)

\[E[f(x)] = \int f(x)\,P(x)\,dx,\]

permette di introdurre i due parametri chiave che quantificano il baricentro e la larghezza di una distribuzione:

• Media definita come \(\mu=E[x]\);
• Varianza definita come \(\sigma^2 = E[(x-\mu)^2]=E[x^2]-\mu^2=\mathrm{Var}(x)\).

Nel caso particolare della distribuzione normale questi due parametri la determinano completamente:

\[P(x) =\frac{\exp\left[-\cfrac{(x-\mu)^2}{2\sigma^2}\right]}{\sqrt{2\pi\sigma^2}}.\]

In questo caso, tramite integrale, si possono quantificare dei ben noti intervalli di confidenza:

intervallo	probabilità	fuori intervallo una volta ogni
\([\mu-\sigma,\ \mu+\sigma]\)	68.27 %	\(\approx 3\) estrazioni
\([\mu-2\sigma,\ \mu+2\sigma]\)	95.45 %	\(\approx 22\) estrazioni
\([\mu-3\sigma,\ \mu+3\sigma]\)	99.73 %	\(\approx 370\) estrazioni

Definizione 2.2: Estensione a più variabili

Se abbiamo un vettore \(\mathbf{x}=(x_1,\cdots, x_n)\) e una distribuzione \(P(\mathbf{x})\), la definizione di aspettazione è estesa a

\[E[f(\mathbf{x})] = \int f(\mathbf{x})\,P(\mathbf{x})\,d^nx,\]

che permette di introdurre

• Vettore medio definito come \(\boldsymbol{\mu}=E[\mathbf{x}]\);
• Matrice di covarianza definita come \(\sigma^2_{ij} = E[(x_i-\mu_i)(x_j-\mu_j)]\).

Due variabili \(x_i\) e \(x_j\) con \(i\neq j\) si dicono indipendenti se la probabilità congiunta è un prodotto \(P(x_i,x_j)=P(x_i)P(x_j)\), e si dicono scorrelate se \(\sigma^2_{ij}=0\). Si può facilmente dimostrare che l’indipendenza implica la scorrelazione; l’implicazione contraria è vera nel caso della gaussiana, ma non in generale.

La naturale estensione della distribuzione normale è la gaussiana multivariata

\[P(\mathbf{x}) = \frac{\exp\left[-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})\cdot\mathcal{P}(\mathbf{x}-\boldsymbol{\mu})\right]}{\sqrt{(2\pi)^n\det\sigma^2}},\]

dove \(\mathcal{P}=\left(\sigma^2\right)^{-1}\) è l’inversa della covarianza, detta matrice di precisione. La naturale estensione dell’intervallo di confidenza \([\mu-\sigma,\mu+\sigma]\), o \((x-\mu)^2/\sigma^2\le1\), è un ellissoide di confidenza definito come

\[(\mathbf{x}-\boldsymbol{\mu})\cdot\mathcal{P}(\mathbf{x}-\boldsymbol{\mu})\le1.\]

Sul significato dei termini della matrice e altri dettagli

In primis, si noti che la probabilità nell’ellissoide non è più circa 68%: il suo valore è inferiore, dipende dal numero di parametri e nel caso di due vale poco meno del 40%. Non lo useremo mai e la prassi sarà la seguente:

• La varianza di \(x_i\) è \(\sigma_i^2=\sigma_{ii}^2\), ossia l’elemento diagonale della matrice di covarianza. Il motivo è il seguente: se consideriamo la probabilità di ottenere \(x_i\), a prescindere dal valore degli altri parametri, dobbiamo integrare le altre variabili in \(P(\mathbf{x})\) e quello che resta è una gaussiana con varianza \(\sigma^2_{ii}\). Graficamente, questo corrisponde a integrare non sull’ellissoide ma su una banda tangente. L’intervallo \([\mu_i-\sigma_i,\mu_i+\sigma_i]\) corrisponde come sempre a circa il 68% di confidenza.
• I termini fuori diagonale rendono l’ellissoide distorto e fanno sì che le fluttuazioni di due variabili \(x_i\) e \(x_j\) non siano casuali ma correlate: se “positivamente” tendono ad avere lo stesso segno, se “negativamente” tendono ad avere segno opposto. Questo può essere quantificato dal coefficiente di correlazione \(\rho=\sigma^2_{ij}/\sigma_i\sigma_j\in[-1,+1]\).

import { renderEllisse } from "./_interactive/ch02-ellisse.js"
{
  const container = document.createElement("div");
  container.style.cssText = "width:100%;position:relative;margin:0 auto;";
  yield container;
  renderEllisse(container);
}

Figura 2.3: L’ellissoide di confidenza su due parametri integra solo il 39.35% della distribuzione, ma se si considera l’incertezza di un parametro a prescindere dal valore degli altri la distribuzione è ancora una Gaussiana standard, con una varianza pari all’elemento diagonale della matrice \(\sigma_{ij}^2\). Il cursore controlla gli elementi fuori diagonale, che distorcono la distribuzione e correlano le fluttuazioni di \(x_1\) e \(x_2\).

TeoremaRegole di propagazione degli errori statistici

Ricordiamo anche le regole di propagazione degli errori statistici: supponiamo che \(y=g(\mathbf{x})\) e di voler determinare, note le fluttuazioni delle singole \(x_i\), quanto fluttui \(y\). Espandendo linearmente \(g\) abbiamo che

\[\sigma_y^2 = \sum_{ij} \frac{\partial g}{\partial x_i}\frac{\partial g}{\partial x_j}\sigma_{ij}^2 = \nabla g \cdot\left( \sigma^2\nabla g\right), \tag{2.1}\]

che, in assenza di correlazioni, diventa la classica somma in quadratura

\[\sigma_y^2 = \sum_i \left(\frac{\partial g}{\partial x_i}\right)^2\sigma_i^2.\]

Dimostrazione

Nell’approssimazione lineare \(g(\mathbf{x}) = g(\boldsymbol{\mu})+\nabla g|_\boldsymbol{\mu}\cdot(\mathbf{x}-\boldsymbol{\mu})\) il valore medio di \(y\) vale

\[\mu_y=E[y] = g(\boldsymbol{\mu})+\nabla g\cdot E[\mathbf{x}-\boldsymbol{\mu}]=g(\boldsymbol{\mu})\]

e quindi, dato che lo scarto vale \(y-\mu_y = \nabla g\cdot(\mathbf{x}-\boldsymbol{\mu})\), la regola sopra citata segue da semplice algebra.

E le derivate superiori? Premesso che se è nota \(P(\mathbf{x})\) si può sempre ipoteticamente fare un integrale di aspettazione, tutto questo formalismo fallisce nella misura in cui le derivate superiori a quella lineare sono importanti. Come esempio banale, si consideri una singola variabile \(x\) a media nulla, \(E[x]=0\), e con una varianza \(\sigma^2\). Se prendiamo \(y=g(x)=x^2\), la regola appena citata predice \(\mu_y=0\) e \(\sigma_y^2=0\). Entrambe le predizioni sono palesemente false: le fluttuazioni esistono e sono tutte sbilanciate verso \(y>0\), quindi la media non è certo zero. Per stimare \(\mu_y\) bisogna andare oltre il limite lineare e calcolare davvero gli integrali, per esempio

\[E[y]=E[x^2]=\sigma^2\neq 0\]

mentre il valore di \(\sigma_y^2=E[y^2]-\mu_y^2=E[x^4]-\sigma^4\) richiede la conoscenza di un momento superiore della distribuzione, \(E[x^4]\), quindi non può essere quantificato (!) senza sapere qualcosa di più sulla distribuzione di \(x\): \(\mu\) e \(\sigma^2\) non bastano a dare una risposta completa.

2.2 Il fit ai minimi quadrati

Ricordiamo le ipotesi che stanno alla base delle procedure di fit. Si suppone di misurare un parametro \(y\) in funzione di un parametro di controllo \(x\), per esempio la corrente in funzione della tensione applicata. La procedura di misura produce una serie di punti sperimentali \((x_k, y_k)\), ognuno con incertezze \(\sigma_k\) sulla base delle caratteristiche di precisione del metodo di misura.

TeoremaMinimizzazione del \(\chi^2\)

Sotto le seguenti ipotesi:

• le misure \(y_k\) sono affette solo da errore statistico di tipo gaussiano;
• l’errore è solo sulle \(y_k\) mentre le \(x_k\) sono note esattamente;
• gli errori dei vari punti sperimentali sono scorrelati;

e dato un ipotetico modello \(f(\mathbf{p}, x)\) con parametri \(\mathbf{p}\), il valore più probabile di \(\mathbf{p}\) è quello che minimizza

\[\chi^2(\mathbf{p}) = \sum_k \frac{\left[y_k - f(\mathbf{p}, x_k)\right]^2}{\sigma_k^2}, \tag{2.2}\]

dove \(w_k = 1/\sigma_k^2\) sono i pesi statistici e \(r_k = y_k - f(\mathbf{p},x_k)\) sono i residui del fit. I parametri \(\mathbf{p}_\mathrm{BF}\) per cui il \(\chi^2\) è minimo sono detti di best fit. Si noti che il best fit non cambia se gli errori vengono moltiplicati per una costante.

Dimostrazione

Date le ipotesi e supponendo che il modello con parametri \(\mathbf{p}\) sia vero, la probabilità di ogni singola acquisizione vale

\[p(x_k,y_k|M) = P_k = \frac{\exp\!\left[-\frac{[y_k - f(\mathbf{p},x_k)]^2}{2\sigma_k^2}\right]}{\sqrt{2\pi\sigma_k^2}},\]

e, dato che gli errori sono indipendenti, la probabilità di ottenere tutta la misura (likelihood) è un prodotto

\[p(D|M) = \prod_k P_k \propto \prod_k \exp\!\left[-\frac{[y_k - f(\mathbf{p},x_k)]^2}{2\sigma_k^2}\right] = e^{-\chi^2/2}, \tag{2.3}\]

dove \(D\) sono i dati \((x_k, y_k)\) e \(M\) il modello \(f(\mathbf{p},x)\). Secondo l’approccio di Bayes, la misura sperimentale aggiorna la probabilità che il modello \(M\) sia vero in base al teorema di Bayes

\[p(M|D) = \frac{p(D|M)\,p(M)}{p(D)} \propto p(D|M) \propto \exp\left(-\frac{\chi^2}{2}\right), \tag{2.4}\]

dove \(p(M)\) è la probabilità “a priori” (prima dell’esperimento, in particolare) che il modello \(M\) potesse essere valido. Assumendo che \(p(M)\) non privilegi alcun modello e notato che \(p(D)\) non dipende da \(M\), minimizzare il \(\chi^2\) equivale a trovare il modello più probabilmente corretto.

2.2.1 L’incertezza del best fit

Ribadiamo: stimare qualcosa senza specificare l’incertezza è sostanzialmente inutile e non andrebbe fatto mai. Sempre sotto le stesse ipotesi alla base del best fit, la distribuzione di probabilità dei parametri di fit è una gaussiana multivariata. Spesso le incertezze vengono prese come l’output di una qualche funzione “black box” di fit, ma non è inutile (né difficile!) rivedere il dettaglio matematico.

TeoremaMatrice di covarianza dei parametri di fit

Definito lo Jacobiano pesato come

\[\mathcal{J}^*_{ij} = \frac{1}{\sigma_i}\left.\frac{\partial f}{\partial p_j}\right|_{x_i}, \tag{2.5}\]

la matrice di covarianza è l’inversa del prodotto fra lo Jacobiano pesato e il suo trasposto

\[\sigma_\mathbf{p}^2 = \left[(\mathcal{J}^*)^T\mathcal{J}^*\right]^{-1}, \tag{2.6}\]

si noti come \(\sigma^2_\mathbf{p}\) sia direttamente proporzionale alla varianza dei dati: se aumenta, \(\sigma_\mathbf{p}^2\) aumenta.

Dimostrazione

La funzione \(\chi^2\) raggiunge un punto di minimo a \(\mathbf{p}_\mathrm{BF}\) quindi il suo gradiente nei parametri è nullo. Facendo uno sviluppo al secondo ordine attorno al punto di best fit, con \(\mathbf{p}=\mathbf{p}_\mathrm{BF}+\Delta \mathbf{p}\) si ottiene

\[\frac{\chi^2(\mathbf{p})}{2} \approx \frac{\chi^2(\mathbf{p}_\mathrm{BF})}{2} + \frac{1}{2}\Delta\mathbf{p}\cdot \mathcal{H}\,\Delta\mathbf{p}, \tag{2.7}\]

dove \(\mathcal{H}\) è l’Hessiano della funzione \(\chi^2(\mathbf{p})/2\). Trascurando (!) le derivate seconde di \(f\), che sono comunque moltiplicate anche per residui \(r_k=y_k-f(\mathbf{p},x_k)\) sperabilmente piccoli, si ottiene

\[ \begin{aligned} \mathcal{H}_{ij} &=\frac{\partial^2}{\partial p_i\,\partial p_j}\frac{\chi^2}{2} = \sum_{k} \frac{1}{\sigma_k^2}\left\{ \left.\frac{\partial f}{\partial p_i}\right|_{x_k} \left.\frac{\partial f}{\partial p_j}\right|_{x_k} - \left.\frac{\partial^2 f}{\partial p_i\,\partial p_j}\right|_{x_k} r_k \right\}\\ &\approx \sum_{k} \frac{1}{\sigma_k^2} \left.\frac{\partial f}{\partial p_i}\right|_{x_k} \left.\frac{\partial f}{\partial p_j}\right|_{x_k} = \sum_k \mathcal{J}^*_{ki}\mathcal{J}^*_{kj} = \left[(\mathcal{J}^*)^T\mathcal{J}^*\right]_{ij} \end{aligned} \tag{2.8}\]

Sotto queste approssimazioni \(p(M|D)\propto\exp(-\chi^2/2)\) diventa una gaussiana multivariata con \(\mathcal{H}\) come matrice di precisione, quindi la matrice di covarianza è l’inversa di \(\mathcal{H}\).

AttenzioneImplementazione in python

Best fit in python. La funzione più generale di best fit disponibile è scipy.optimize.curve_fit, che restituisce: i parametri di best fit e la matrice di covarianza \(\sigma_p^2\). La sintassi è la seguente

# Esempio di codice per fit con una funzione p1*sin(x)+p2
#
# [Argomenti]
# xv    vettore x
# yv    vettore y
# p0    parametri iniziali
# s0    vettore degli errori
#
# [Risultati]
# pp    best fit
# CovB  matrice di covarianza
p0 = [1, 0]
pp, CovB = curve_fit(lambda x, p1, p2: p1*np.sin(x) + p2, xv, yv, p0=p0, sigma=s0, absolute_sigma=True)
sigma_p = np.sqrt(np.diag(CovB))        

e gli errori vengono intesi come assoluti, ossia l’equivalente di absolute_sigma = True. Se si omette il parametro absolute_sigma o lo si imposta a False, che è il comportamento standard della funzione, gli errori vengono invece riscalati come descritto nel riquadro seguente.

AttenzioneImplementazione in MATLAB e absolute_sigma

Best fit in MATLAB. La funzione più generale di best fit disponibile è nlinfit, che restituisce: i parametri di best fit, lo Jacobiano, i residui, la matrice di covarianza \(\sigma_p^2\) e il Mean Squared Error (MSE), che non è altro che il \(\chi^2\) ridotto, \(\chi^2_\nu\). La funzione nlinfit non implementa alcun equivalente del parametro absolute_sigma e riscala sempre (!) gli errori. La sintassi è la seguente

% Esempio di codice per fit con una funzione p1*sin(x)+p2
%
% [Argomenti]
% xv    vettore x
% yv    vettore y
% fitf  funzione di fit (inline qui)
% beta0 parametri iniziali
% wv    vettore dei pesi 1/sigma^2
%
% [Risultati]
% beta  best fit
% R     vettore residui
% J     matrice Jacobiana
% CovB  matrice di covarianza
% MSE   Mean Square Error, o anche chi quadro ridotto 

beta0 = [1 0];
[beta, R, J, CovB, MSE] = nlinfit(xv, yv, @(p,x) p(1)*sin(x) + p(2), beta0, Weights = wv);
sigma_p = sqrt(diag(CovB));         % equivalente di absolute_sigma = False

Implementare l’equivalente di absolute_sigma = True è tuttavia banale e l’operazione di riscalatura degli errori è facilmente invertibile. Gli errori vengono moltiplicati per un fattore \(\alpha\) che rende il \(\chi^2_\nu\) uguale a 1

\[\begin{aligned} \sigma_k &\to \alpha\sigma_k\\ \chi^2_\nu &\to \chi^2_\nu/\alpha^2 = 1\\ \sigma^2_\mathbf{p} &\to \sigma^2_\mathbf{p}|_{res}=\alpha^2\sigma^2_\mathbf{p} \end{aligned} \tag{2.9}\]

da cui la riscalatura della covarianza va fatta con \(\alpha^2=\chi^2_\nu\): per invertirla e tornare agli errori assoluti basta dividere la matrice di covarianza per \(\chi^2_\nu\). Quindi rispetto allo snippet precedente dobbiamo calcolare gli errori come segue

sigma_p = sqrt(diag(CovB / MSE));   % equivalente di absolute_sigma = True

Vero o falso. In generale vogliamo che, a meno che per qualche motivo sia impossibile quantificare gli errori statistici dei dati, usiate l’equivalente di absolute_sigma = True e prendiate le incertezze \(\sigma_k\) per quelle che sono, senza ammettere riscalatura. Questo è fondamentale per poter fare test statistici chiave sul fit, come quello dei residui o del \(\chi^2_\nu\). Nel caso in cui il test fallisca, daremo per scontato che le incertezze che ne derivano possano essere completamente inaffidabili (più dettagli in seguito).

Errori sistematici dei dati. Ricordiamo che questi non vanno inclusi nelle \(\sigma_i\), perché rompono le ipotesi di fit in vario modo: non sono statistici, non sono gaussiani, sono completamente correlati. Per chiarire poi quanto la cosa sia problematica si consideri un semplice partitore resistivo che soddisfa la legge

\[V_{out} = \frac{R_1}{R_1+R_2}V_{in},\]

in cui ci immaginiamo di variare \(V_{in}\), misurando sia \(V_{in}\) che \(V_{out}\) e quindi fittando la slope della curva \(V_{out}(V_{in})\) per stimare \(R_2/R_1\). Il misuratore di voltaggio avrà sicuramente dei limiti di accuratezza, ma da cosa derivano? Un errore sistematico di offset è chiaramente irrilevante dato che non impatta la slope. In maniera simile un errore di scala può avere un impatto solo nella misura in cui è diverso fra il canale che misura \(V_{out}\) e quello che misura \(V_{in}\). In conclusione, includere l’accuratezza dei canali di misura gonfierebbe in maniera ingiustificata le incertezze. Quello che va fatto è fittare i dati con i soli errori statistici e solo dopo chiedersi come una accuratezza finita possa indurre l’errore di stima che ne deriva. Come spesso succede, per gli errori sistematici non ci sarà una ricetta precisa sempre applicabile.

AttenzioneRegressione lineare standard

Le formule date per la matrice di covarianza \(\sigma^2_\mathbf{p}\) si applicano perfettamente al caso della regressione lineare, ovviamente, dato che la linearizzazione è un dato di fatto e non una approssimazione. Qui per non appesantire troppo assumiamo che gli errori dei dati siano tutti uguali, quindi fisseremo \(\sigma_k=\sigma\). Data la funzione di fit

\[f(m,q,x) = mx+q,\]

la matrice Jacobiana è definita come

\[\mathcal{J}_{ij} = \left.\frac{\partial f}{\partial p_j}\right|_{x_i} = \begin{pmatrix} x_1 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{pmatrix}\]

mentre lo Jacobiano pesato è \(\mathcal{J}^*=\mathcal{J}/\sigma\). Grazie alla linearità in questo caso abbiamo la gradevole relazione esatta \(\left\{f(\mathbf{p},x_1),\cdots,f(\mathbf{p},x_n)\right\} = \mathcal{J}\mathbf{p}\), dove \(\mathbf{p}=(m,q)\) è il vettore dei parametri.

Best fit. Partendo dalla definizione di \(\chi^2\) e chiamando \(\mathbf{y}=(y_1,\cdots,y_n)\) il vettore delle misure

\[\chi^2(\mathbf{p}) = \frac{|\mathbf{y}-\mathcal{J}\mathbf{p}|^2}{\sigma^2}\]

con un punto stazionario che è banale da determinare (è una parabola in due dimensioni) e richiede

\[\mathcal{J}^T(\mathbf{y}-\mathcal{J}\mathbf{p}_\mathrm{BF}) = \mathbf{0},\]

da cui

\[\mathbf{p}_\mathrm{BF} = \left(\mathcal{J}^T\mathcal{J}\right)^{-1}\mathcal{J}^T\mathbf{y}.\]

Questa formula può essere facilmente esplicitata, considerando la matrice \(\mathcal{J}^T\mathcal{J}\) e la sua inversa

\[ \begin{aligned} \mathcal{J}^T\mathcal{J} &= \begin{pmatrix} \sum x^2_k & \sum x_k \\ \sum x_k & n \end{pmatrix}= \begin{pmatrix} S_{xx} & S_x \\ S_x & n \end{pmatrix}\\ \left(\mathcal{J}^T\mathcal{J}\right)^{-1} &= \frac{1}{\Delta} \begin{pmatrix} n & -S_x \\ -S_x & S_{xx} \end{pmatrix} \end{aligned} \]

dove \(S_x=\sum x_k\) e \(S_{xx}=\sum x_k^2\) e il determinante \(\Delta = nS_{xx}-S_x^2\). A questo punto

\[ \mathbf{p}_\mathrm{BF} = \frac{1}{\Delta} \begin{pmatrix} n & -S_x \\ -S_x & S_{xx} \end{pmatrix} \mathcal{J}^T\mathbf{y} = \frac{1}{\Delta} \begin{pmatrix} n & -S_x \\ -S_x & S_{xx} \end{pmatrix} \begin{pmatrix} S_{xy} \\ S_y \end{pmatrix} = \frac{1}{\Delta} \begin{pmatrix} nS_{xy}-S_xS_y \\ S_yS_{xx}-S_xS_{xy} \end{pmatrix}. \]

Covarianza di fit. La matrice di precisione vale (stesso prodotto già visto a parte un peso \(1/\sigma^2\))

\[ \mathcal{P} = \left(\mathcal{J}^*\right)^T\mathcal{J}^* = \frac{1}{\sigma^2} \begin{pmatrix} S_{xx} & S_x \\ S_x & n \end{pmatrix}, \]

e infine le incertezze di fit sono descritte da

\[ \sigma^2_\mathbf{p} = \mathcal{P}^{-1} =\frac{\sigma^2}{\Delta}\begin{pmatrix} n & -S_x \\ -S_x & S_{xx} \end{pmatrix}. \]

Si noti che la presenza di una correlazione positiva o negativa fra le incertezze di \(m\) e \(q\) dipende dal segno di \(\sum x_k\) ossia dal baricentro dei punti fittati: se è localizzato ad \(x\) positivi allora la correlazione è negativa e viceversa.

Come esercizio opzionale, si noti che una ridefinizione dell’origine delle \(x\) non cambia il valore di \(\Delta\) e quindi neanche di \(\sigma_m^2\). Impatta invece sia su \(S_x\) che su \(S_{xx}\), e quindi sull’incertezza \(\sigma^2_q\) e sulla correlazione \(\sigma^2_{mq}\). Tutto questo è consistente con il fatto che una traslazione dell’origine delle \(x\) di un dato \(\Delta x\) cambia la definizione di intercetta \(q\to q+m\Delta x\), che quindi può “ereditare” una incertezza dal parametro \(m\) in base alle regole di propagazione degli errori.

2.2.2 Fit playground

Ricapitoliamo qui alcune problematiche connesse con le procedure di fit con il supporto di Figura 2.4, che può generare una sequenza di fit casuali, per visualizzare alcuni concetti e fenomenologie di base in un fit.

Significato della matrice di covarianza. La matrice \(\sigma_\mathbf{p}^2\), se calcolata a partire dai veri errori \(\sigma_i\), quantifica l’incertezza statistica del fit. Se pensiamo ad acquisizione e fit come una misura di \(\mathbf{p}\), allora \(\sigma_\mathbf{p}^2\) è semplicemente l’incertezza statistica di questa misura, propagata dai dati. Dato che le incertezze statistiche si possono quantificare sperimentalmente, oltre che propagare, Figura 2.4 permette di verificare che ripetendo il fit molte volte e stimando dai best fit stessi la covarianza, si ottengono effettivamente incertezze coincidenti.

Covarianza e correlazione. I parametri \(m\) e \(q\) del fit lineare possono o meno essere correlati: in un fit lineare ad errori sui dati costanti, come in questo caso, questo dipende solo dal baricentro \(x\) dei dati. Nella Figura 2.4 è possibile trascinare i punti portandoli verso \(x>0\) o \(x<0\), e apprezzare come l’ellissoide di confidenza si allunghi in una direzione. Il significato della correlazione è ovvio: se per esempio i dati sono a destra dell’origine \(x\), ogni fluttuazione della stima di \(q\) richiede una correzione di segno opposto di \(m\) per avere un fit che passi comunque per i dati sperimentali.

import { renderFitCovariance } from "./_interactive/ch02-fit-covarianza.js"
{
  const container = document.createElement("div");
  container.style.cssText = "width:100%;position:relative;margin:0 auto;";
  yield container;
  renderFitCovariance(container);
}

Figura 2.4: Simulatore di fit. Sinistra: dataset sintetico in base ad un modello lineare più rumore gaussiano, gli slider permettono anche di inserire un salto fra il ramo positivo e negativo non incluso nel modello, che simula una possibile non-idealità non prevista del sistema di misura; in basso i residui normalizzati. Destra: parametri di fit con ellisse di confidenza a una \(\sigma\) nel piano \((m,q)\); il grafico include anche una stima dell’ellisse di confidenza sulla base dei risultati di fit sintetici ripetuti. In basso è riportato un istogramma dei \(\chi^2_\nu\) ottenuti, confrontato con la distribuzione attesa.

Quando il fit fallisce. Esploriamo ora aspetti più problematici e fit non ideali, partendo dall’assunto che vogliamo in ogni caso partire da dei fit fatti con gli errori statistici veri dei dati. Come noto, la qualità di fit può essere quantificata con il \(\chi^2_\nu = \chi^2 / \nu\) con \(\nu = N - P\), dove \(N\) è il numero di dati e \(P\) il numero di parametri di fit. Se il modello è corretto, il valore atteso è circa 1 con una distribuzione specifica, la cui larghezza per \(\nu\) grande tende a \(\sqrt{2/\nu}\).

Qualsiasi strumento di misura se spinto al limite mostra delle non-idealità e nella misura sintetica di Figura 2.4 questo è implementato con un “salto” \(\Delta\) aggiunto fra le misure per \(x\) positive e negative. Immaginandoci che lo sperimentatore non sia a conoscenza di questo problema, il fit non contiene questa correzione e quindi stima “male” la pendenza della retta: aumentando \(\Delta\) il best fit si sposta sistematicamente verso \(m\) inferiori alla pendenza vera. Una osservazione chiave è che \(\sigma_\mathbf{p}^2\) continua a descrivere correttamente la fluttuazione statistica dei fit, ma non traccia in nessun modo l’emergere di questo problema di stima. Si noti che anche arrivando a \(\Delta=6\) con \(\sigma=2\) la deviazione del \(\chi^2_\nu\) non è nemmeno ovvia da notare.

Calando la \(\sigma\) della misura sintetica (nel mondo reale questo potrebbe succedere facendo molte medie o prendendo molti punti!) il \(\chi^2_\nu\) aumenta, finché la presenza di un problema diventa inequivocabile, per esempio arrivando a \(\sigma=0.2\). Due osservazioni cruciali qui:

• a prescindere da quanto sia grande \(\sigma\) e a prescindere da quanto sia brutto \(\chi^2_\nu\) la \(\sigma_\mathbf{p}^2\) continua a fare il suo lavoro, e stima correttamente la fluttuazione statistica della procedura di best fit;
• quando \(\sigma\) diventa molto piccolo, \(\sigma_\mathbf{p}^2\) diventa una stima insensata dell’incertezza: il fit è effettivamente sempre più preciso, ma manca malamente il valore vero; non è sempre così, ma questo esempio mostra che può tranquillamente succedere.

A questo punto esistono le seguenti opzioni:

Opzione 0: reality check e ragionevolezza

Se la problematica è emersa perché si è deciso di fare un numero irragionevolmente grande di medie, acquisendo un numero irragionevolmente grande di punti sperimentali, magari per esempio cercando di studiare la risposta di un filtro passa banda audio campionandolo da \(10\,\mathrm{mHz}\) fino a \(10\,\mathrm{MHz}\) e poi pretendendo di fittare tutto perfettamente e in un colpo solo… ci si dovrebbe seriamente chiedere quale sia l’obiettivo. In questo ricordiamo che se degli effetti sistematici mantengono alcuni residui a valori finiti, l’impatto sul \(\chi^2_\nu\) è proporzionale a \(1/\sigma^2\), quindi può solo peggiorare mediando di più. Inoltre, la scala di deviazione tollerabile per \(\chi^2_\nu\) è \(\sqrt{2/\nu}\) e diventa sempre più stringente più aumentano i punti di misura. Quindi a meno che il vostro obiettivo non sia esplicitamente di scovare tutti i possibili limiti del metodo sperimentale adottato (ed esistono, sempre), qualcosa non torna nelle scelte fatte.

Ergo, se questo è il caso, il consiglio è di rifare la misura in condizioni sperimentali più sagge.

Opzione 1: risolvere il problema

Un \(\chi^2_\nu\) molto fuori range non è necessariamente un problema, è un segnale: un segnale che qualcosa non va come ci aspettiamo. Questo è magari perfino interessante perché magari ci porta a capire qualcosa di nuovo. Un primo intento quindi dovrebbe essere di osservare i residui, e cercare di capire se c’è un problema di modello, o magari un problema di misura. O forse anche un problema di stima delle incertezze dei dati, se per esempio i residui non sono correlati. Come unico caveat sulla modifica del modello di fit, è chiaro che aggiungendo infiniti parametri si può fittare qualsiasi cosa, ma questo di nuovo è un non-obiettivo: se si vuole avere una curva che passa perfettamente per i punti sperimentali una buona opzione è sempre quella di prendere una penna e unire i punti come sulla Settimana Enigmistica, mentre qui l’obiettivo dovrebbe essere diverso e vorremmo invece capire/stimare qualcosa.

Opzione 2: fare quel che si può

Esistono varie strategie. In primis, collegandosi al caso specifico del filtro, non è assolutamente necessario fittare la banda passante fino a soppressioni insensatamente piccole e su range spettrali insensatamente ampi, estendendo così tanto lo span di frequenze da rendere probabilmente inadeguato il modello di fit. Una possibilità è quindi fare un fit su un dominio più ridotto in cui il modello funziona e campiona bene il parametro target, magari verificando se i risultati sono stabili con la scelta del dominio o emergono deviazioni sistematiche (ora sì!) quantificabili. Restando per esempio sullo specifico esempio del filtro, dato che l’obiettivo è di determinare i bordi del filtro, si può perfino semplicemente fittare i singoli roll-off. Rimandiamo alle Note di Lab1 per ulteriori idee.

Anche se non è certo il primo consiglio e va fatto come extrema ratio dopo aver tentato tutto il resto, verrà tollerato l’utilizzo di \(\sigma^2_\mathbf{p}|_{res}\), ossia il risultato nell’equivalente della configurazione absolute_sigma = False. Deve tuttavia essere chiaro quali siano il senso, i limiti e le fallacità di questa scelta. Nella Figura 2.4 \(\sigma^2_\mathbf{p}|_{res}\) è mostrato con l’ellissoide tratteggiato: a \(\Delta=6\) si vede come allo scendere di \(\sigma\) quantomeno ha pragmaticamente il merito di non diventare insensatamente piccolo come \(\sigma^2_\mathbf{p}\); il suo valore sostanzialmente satura alla scala di incertezza statistica che era presente per la \(\sigma\) a cui \(\chi^2_\nu\) ha iniziato a deviare significativamente da \(1\), e quindi a cui è diventato chiaro che il fit fosse problematico. In questo senso, è una scelta conservativa che ha un suo senso. Resta tuttavia chiaro dall’esempio in Figura 2.4 che la discrepanza dal valore vero è comunque molte volte più grande anche di questa incertezza, che infatti non ha una chiara valenza statistica.

Si noti che un approccio simile, anche se in un contesto non del tutto riconducibile a questo, viene usato in metrologia e va sotto il nome di rapporto di Birge (Birge 1932), che viene attualmente usato per quantificare l’incertezza della costante di gravitazione \(G\) sulla base di diverse fonti sperimentali non del tutto consistenti (Mohr et al. 2025; Bodnar e Elster 2014). Le falle nell’utilizzo di tutto questo dentro una procedura di fit con un \(\chi^2_\nu\) elevato, in particolare in presenza di un modello fallace, sono palesi: la più grave segue dalla forte correlazione dei residui che indica errori non indipendenti, quindi il \(\chi^2\) non ha una forma matematicamente corretta, dato che deriva da un prodotto di probabilità, che ha senso solo per probabilità indipendenti. Quindi il problema diventa sostanzialmente non ben posto.

Come messaggio finale, la cosa più fondamentale di tutte è documentare le scelte fatte e le procedure seguite. Documentare serve anche a permettere a chi viene dopo di approcciarsi al vostro lavoro con la speranza di capirci qualcosa. I destinatari potreste essere voi stessi qualche mese nel futuro, quando probabilmente avrete rimosso dalla memoria la gran parte dei dettagli.

2.2.3 Fit non lineari e altri problemi

Minimi multipli e condizioni iniziali. Nella regressione lineare \(\chi^2(\mathbf{p})\) è sostanzialmente una parabola in più dimensioni, con un unico minimo facilmente individuabile. Questo non è vero quando \(f(\mathbf{p},x)\) è una funzione di fit non lineare. In questo caso \(\chi^2(\mathbf{p})\) può avere molti minimi relativi, su cui la procedura di fit può facilmente “incagliarsi”. Per questo motivo è fondamentale fare dei test con i parametri di fit e partire da una condizione sufficientemente “vicina” al best fit cercato. La Figura 2.5 illustra il problema nel caso di un segnale periodico \(y=\sin(2\pi ft+\varphi)\) e mostra come la traiettoria di minimizzazione possa portare a diversi risultati in base al guess iniziale per i parametri.

import { renderFitNonlineare } from "./_interactive/ch02-fit-nonlineare.js"
{
  const container = document.createElement("div");
  container.style.cssText = "width:100%;position:relative;margin:0 auto;";
  yield container;
  renderFitNonlineare(container);
}

Figura 2.5: Esempio di fit non lineare a due parametri. Sinistra: una sinusoide di ampiezza nota campionata con pochi punti per periodo; in verde il modello coi parametri selezionati, in rosso tratteggiato la legge vera. In basso i residui normalizzati. Il parametro \(\Delta\) introduce un offset fra i dati a \(x\) positiva e negativa, simulando una imperfezione sperimentale ignota. Destra: la funzione \(\chi^2_\nu(f, \varphi)\) ha vari minimi, alcuni dovuti alla periodicità della fase, altri al fatto che per via dell’aliasing lo stesso dataset può essere fittato in maniera ottimale usando due frequenze diverse.

Messaggi di fit misteriosi: Jacobiano quasi singolare. Chiudiamo con un commento su una problematica su cui presto o tardi capita di scontrarsi: un fit che non converge. La non convergenza di un fit è assicurata se il minimo del \(\chi^2\) non è un punto, bensì una retta o una curva: in questo caso la procedura di fit non può dare una risposta univoca al best fit. Per ottenere questa spiacevole situazione è sufficiente inserire un parametro di fit che non ha nessun ruolo. In questo caso diciamo che lo Jacobiano è singolare e infatti una delle colonne di \(\mathcal{J}\) è completamente nulla.

Una situazione non molto diversa si ottiene anche quando la funzione di fit \(f(\mathbf{p},x)\) ha semplicemente una dipendenza molto debole da uno dei parametri, che chiamiamo \(p\). In questo caso sebbene lo Jacobiano non sia esattamente singolare, la procedura numerica di minimizzazione del \(\chi^2\) può faticare molto a convergere. Sebbene il tutto dipenda molto dalla funzione di fit specifica e dall’algoritmo di minimizzazione del \(\chi^2\), il consiglio per non incappare in problemi numerici è di evitare di fittare con parametri che differiscono per vari ordini di grandezza nello stesso fit. Per esempio, se la risposta di un filtro è controllata da un tempo di rilassamento di \(\tau\approx10^{-2}\,\mathrm{s}\), con magari una risonanza attesa a \(\nu_0\approx 10^6\,\mathrm{Hz}\), usare \(\tau\) e \(\omega_0\) come parametri di fit potrebbe essere incauto dato che i due best fit attesi differiscono di 8 ordini di grandezza. In questo caso basta fittare in \(\gamma=1/\tau\approx 10^2\) o magari cambiare unità per \(\nu_0\) ed esprimerlo in MHz.

In conclusione…

ImportanteTake-home

Alla fine di questo recap è fondamentale avere ripreso confidenza con i seguenti concetti:

• Cosa sono precisione e accuratezza
• Cosa è una matrice di covarianza e come è connessa con l’estensione a più parametri della distribuzione normale.
• Regole di propagazione degli errori statistici, inclusa covarianza.
• Procedura di best fit e aspetti subdoli: quali sono le ipotesi, in che modo un fit non lineare è peculiare? qual è il significato dell’incertezza di fit?
• Gestione dei fit patologici e comprensione della differenza fra usare l’errore dei dati in senso assoluto o relativo.

Birge, Raymond T. 1932. «The Calculation of Errors by the Method of Least Squares». Physical Review 40: 207–27. https://doi.org/10.1103/PhysRev.40.207.

Bodnar, Olha, e Clemens Elster. 2014. «On the Adjustment of Inconsistent Data Using the Birge Ratio». Metrologia 51 (5): 516–21. https://doi.org/10.1088/0026-1394/51/5/516.

JCGM. 2012. International Vocabulary of Metrology — Basic and General Concepts and Associated Terms (VIM). JCGM 200:2012. 3ª ed. Joint Committee for Guides in Metrology. https://doi.org/10.59161/JCGM200-2012.

Mohr, Peter J., David B. Newell, Barry N. Taylor, e Eite Tiesinga. 2025. «CODATA Recommended Values of the Fundamental Physical Constants: 2022». Reviews of Modern Physics 97 (2): 025002. https://doi.org/10.1103/RevModPhys.97.025002.