Esperienze di fisica

L'analisi di regressione e il metodo dei minimi quadrati

Introduzione

Il problema della ricerca di una relazione funzionale tra due o più insiemi di dati in modo che si possa stabilire tra essi una qualche relazione matematica costituisce un aspetto standard nell'ambito delle scienze fisiche e statistiche. Con l'odierna grande disponibilità di dati provenienti da ambiti diversi quali le scienze biologiche, l'economia, le scienze sociali, il management, una tale esigenza si è indubbiamente estesa e si è fatta più urgente la necessità di sintetizzare la massa dei dati estraendone, per mezzo di tecniche più o meno automatizzate (oggi individuate dal termine data mining), relazioni tra le grandezze che li caratterizzano. L'analisi di regressione o semplicemente, regressione, utilizza le tipiche tecniche della statistica e risponde a questa esigenza fornendo validi strumenti per modellizzare le relazioni esistenti tra grandezze e, probabilmente, rappresenta la tecnica statistica più utilizzata.

Con il supporto di semplici visualizzazioni intendiamo sviluppare in questa pagina (regressione01.html) un percorso didattico di carattere teorico che presenti gli aspetti di base del metodo dei minimi quadrati tramite il quale si giunge ad individuare quella particolare relazione in grado di adattarsi nel modo migliore ad un dato insieme di osservazioni sperimentali. Nella pagina che segue, regressione02.html, presentiamo invece esempi della applicazione di queste tecniche in un laboratorio didattico.

Dividiamo quanto segue in due parti in quanto nella prima trattiamo il metodo dei minimi quadrati con un costante riferimento ad una funzione lineare mentre nella seconda parte estendiamo il metodo a funzioni più generali.

Poiché ci si riferisce spesso a funzioni di più variabili, come requisito di massima supponiamo siano note le nozioni di analisi necessarie per la ricerca di massimi e/o minimi di funzioni di una variabile e lasciamo invece a livello intuitivo la loro naturale estensione al caso di più variabili. Le dimostrazioni più impegnative sono inoltre riportate in riquadro.

PARTE 1

PARTE 2

Bibliografia

PARTE 1

Proprietà della media aritmetica

In questa sezione ma pure nelle seguenti, supponiamo di disporre di un insieme \(A\) di \(n\) elementi ciascuno dei quali è costituito da una coppia \((x_i,y_i)\) \begin{equation} A=\Bigl\{(x_1,y_1),(x_2,y_2),\dots (x_{n-1},y_{n-1}),(x_n,y_n)\Bigr\} \tag{1}\label{eq:01} \end{equation} dove i valori \(x_i\) sono rappresentativi delle misure sperimentali di una grandezza \(X\), mentre i valori \(y_i\) sono le corrispondenti misure della grandezza \(Y\). A prescindere dalle singole misure e quindi in generale, intenderemo che al verificarsi comunque della causa \(x\) (la variabile indipendente) sia \(y\) (la variabile dipendente) la risposta o l'effetto di un certo sistema fisico.

L'insieme di misure \(A\) si rappresenta graficamente tramite un grafico a dispersione o scatter plot (fig. 1) dove i valori della variabile indipendente sono disposti in ascissa e quelli della variabile dipendente in ordinata. Per l'origine sperimentale di queste coppie, ciascun valore è comunque soggetto ad incertezze ed errori ma inizialmente supporremo che tali incertezze possano essere trascurate rispetto ai valori misurati o al loro numero: considereremo pertanto queste coppie prive di errori cosicché nel diagramma a dispersione potremo rappresentarle come semplici punti.

Volendo sintetizzare i valori assunti dalla grandezza \(Y\) con un qualche indice, diverse sono le possibili scelte: potremmo per esempio calcolare la media aritmetica dei valori rilevati \begin{equation} \overline{y}={1\over n}{\sum_{i\,=\,1}^ny_i} \tag{2}\label{eq:02} \end{equation} oppure, in alternativa, potremmo determinare l'intervallo che comprende tutte le misure \(y_i\) e la sua ampiezza \[ [y_{min},y_{max}],\qquad \Delta y=y_{max}-y_{min} \] essendo \(y_{min}\) e \(y_{max}\) il minimo e il massimo dei valori assunti dalla \(Y\). E come ulteriore possibilità potremmo pure scegliere il valore medio tra gli estremi \[ y_{med}= {y_{max}+y_{min}\over 2}. \]

Tra queste grandezze viene comunque preferita la media aritmetica in quanto possiede la notevole proprietà che sotto dimostriamo e che ritroveremo anche in altri contesti. Difatti, definito come scarto dalla media di una misura, la differenza tra la misura stessa e il valore medio \(\overline{y}\) \begin{equation} s_i=y_i-\overline{y}, \tag{3}\label{eq:03} \end{equation} l'insieme degli scarti è costituto sia da valori positivi che eccedono la media, sia da valori negativi. Considerando la loro somma \begin{equation} \sum_{i\,=\,1}^n s_i=\sum_{i\,=\,1}^n(y_i-\overline{y}) \end{equation} e osservato che \(\sum_{i=1}^n\overline{y}=n\,\overline{y}\), riscriviamo la precedente come \[ \eqalign{\sum_{i\,=\,1}^n s_i &=\sum_{i\,=\,1}^n y_i- \sum_{i\,=\,1}^n \overline{y}\cr & =\sum_{i\,=\,1}^n y_i-n\,\overline{y}.} \] In base alla definizione \eqref{eq:02}, non ci resta che sostituire la sommatoria \(\sum_{i\,=\,1}^n \overline{y}\) con \(n\overline{y}\) e troviamo \begin{equation} \sum_{i\,=\,1}^n s_i=n\overline{y}-n\overline{y}=0\colon\tag{4}\label{eq:04} \end{equation} la media aritmetica rappresenta quindi quella grandezza che assegna la medesima importanza agli scarti positivi e negativi e cioè alle misure maggiori o minori del suo valore. Nella figura 1, sovrapponiamo alla distribuzione delle misure la loro media e rappresentiamo gli scarti positivi e negativi come segmenti di diverso colore: la loro somma è comunque nulla.

Fig. 1. Media di un insieme di misure e scarti dalla media.

Modello lineare, minimi quadrati e regressione

Tralasciando per ora il problema di assegnare alla media una qualche stima della sua incertezza e, fissandoci invece sulla ricerca di una relazione che descriva le misure dell'insieme \eqref{eq:01}, osserviamo come in figura 1 ad un andamento crescente della variabile \(x\) si possa far corrispondere un andamento generalmente crescente della variabile \(y\). La retta orizzontale associata al valor medio non sembra quindi adatta a rappresentare adeguatamente questi dati. Il modello che invece intendiamo seguire e che sembra più aderente all'insieme dei dati, si fonda sulla funzione matematica \begin{equation} y=f(x)=a x+b. \tag{5}\label{eq:05} \end{equation} che costituisce una semplice equazione di primo grado e rappresenta la più elementare relazione tra le variabili \(x\) e \(y\). Tale equazione esprime pure una funzione lineare e il suo grafico cartesiano è nient'altro che una linea retta. Per tale motivo il modello scelto per descrivere l'insieme \(A\) verrà indicato come il modello lineare.
È pure evidente che vi possono essere molte altre scelte alternative alla funzione \eqref{eq:05} per cui, come ribadiremo più avanti, il modello che viene scelto in sostituzione della \eqref{eq:05} dovrà essere generalmente formulato sulla base di motivate considerazioni fisiche di carattere teorico o sperimentale.

In corrispondenza della generica coppia \((x_i,y_i)\) possiamo associare il punto che possiede la medesima ascissa ma che giace sulla retta \eqref{eq:05}. Questo è caratterizzato dalle coordinate \((x_i,a x_i+b)\) e potremmo considerare la sua ordinata \(a x_i+b\) come la "previsione" o il "valore atteso" del modello dovuti al verificarsi di \(x_i\). Definita come "deviazione" o "discrepanza" o anche "residuo" o "scarto" del modello la distanza con segno \(e_i\) di questi due punti (fig. 2) \begin{equation} \eqalign{e_i&=y_i-(a x_i+b)\cr &=y_i-a x_i-b\qquad i=1,2,\dots,n,\cr}\tag{6}\label{eq:06} \end{equation}

Fig. 2. Retta di regressione e scarti.

evidentemente questa grandezza avrà segno positivo per quei punti che giacciono nel semipiano superiore della retta (per ora incognita) e segno negativo se al di sotto di essa o nullo se la coppia in questione appartiene alla retta stessa. Volendo comunque introdurre una funzione che possa, in qualche modo, esprimere una qualche distanza della retta dall'intero insieme di osservazioni, potremmo estendere la condizione già applicata agli scarti dalla media, ossia imporre che la funzione di due variabili \(s(a,b\)) definita dalla somma \begin{equation} s(a,b)=\sum_{i\,=\,1}^n e_i=\sum_{i\,=\,1}^n(y_i-ax_i-b).\tag{7}\label{eq:07} \end{equation} abbia valore nullo ossia sia soddisfatta l'equazione \begin{equation} s(a,b)=0.\tag{8}\label{eq:08} \end{equation}

Una tale scelta comunque comporta due problemi il primo dei quali è di carattere algebrico: essendo \(a\) e \(b\) parametri incogniti, la condizione \eqref{eq:08} non permette di individuarli entrambi per cui è necessario aggiungere una qualche altra condizione. Il secondo problema consiste nel fatto che la stessa condizione non identifica univocamente la retta ottimale in quanto possono esistere rette che la soddisfano ma che non descrivono adeguatamente l'insieme delle misure come esemplificato nella figura 3.

Fig. 3. Somma degli scarti nulla ma retta inadeguata.

Scegliamo quindi di sostituire la condizione \eqref{eq:08} con la somma dei quadrati delle differenze \eqref{eq:06}. In tal modo la funzione che definisce la "distanza" della retta dall'insieme delle misure assume la forma \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ s^2(a,b)=\sum_{i\,=\,1}^n e_i^2=\sum_{i\,=\, 1}^n\,(y_i-ax_i-b)^2.\tag{9}\label{eq:09} } \end{equation} Questa espressione, identificata con vari nomi quali funzione obiettivo o, in contesti economici, come funzione costo, dipende come evidenziato a primo membro, dai parametri \(a\) e \(b\) e possiede segno positivo o nullo. Per determinare i valori ottimali \(\hat{a}\) e \(\hat{b}\) dei parametri e di conseguenza individuare quella particolare retta che meglio si adatta ai dati sperimentali, seguiremo il seguente criterio:

in corrispondenza dei valori ottimali dei parametri, \(\hat{a}\) e \(\hat{b}\), la funzione obiettivo presenta il minimo assoluto.

Comprendiamo così come tale criterio applicato alla \eqref{eq:09} dia ragione del nome per cui è noto in ambito scientifico: il metodo dei minimi quadrati consiste pertanto nel minimizzare la funzione obiettivo, somma dei quadrati delle differenze tra osservazioni e modello, così da determinare le grandezze da cui essa dipende e che compaiono come parametri nella sua espressione analitica. L'applicazione di questo criterio nel caso specifico della funzione \eqref{eq:09} identifica la regressione lineare o linear fit mentre la funzione che sta alla base del modello, nel nostro caso la \eqref{eq:05}, è la cosiddetta funzione modello o funzione di fit.

Le condizioni matematiche per individuare il minimo della funzione di due variabili \(s^2(a,b)\) sono una naturale estensione delle tecniche applicabili a funzioni di una variabile e, in questo caso, consistono nell'annullare contemporaneamente le due derivate prime \begin{equation} \cases{\displaystyle {\partial(s^2)\over \partial a}=0\cr\\ \displaystyle{\partial(s^2)\over \partial b}=0\cr} \end{equation} dove si è introdotto per la prima volta il simbolo \(\partial\) e l'operazione di derivata parziale. Supposta quindi una funzione \(z(y_1,y_2,\dots y_n)\) dipendente dalle variabili \(y_1,y_2,\dots,y_n\), la scrittura \(\partial z/\partial y_i\) rappresenta la derivata prima parziale della \(z\) eseguita considerando la sola \(y_i\) come grandezza variabile mentre le rimanenti grandezze \(y_j\) con \(j\not=i\) vanno supposte costanti. Eseguendo il calcolo delle precedenti due derivate otteniamo il sistema di equazioni (da questo punto in poi omettiamo gli indici di sommatoria sostituendo \(\sum_{i\,=\,1}^n\) con \(\sum\)) \begin{equation} \cases{\displaystyle -2\sum x_i(y_i-ax_i-b)=0\cr\\ \displaystyle-2\sum(y_i-a x_i-b)=0\cr}\qquad\Longrightarrow\qquad \cases{\displaystyle \sum x_i(y_i-ax_i-b)=0\cr\\ \displaystyle\sum (y_i-a x_i-b)=0.\cr} \tag{10}\label{eq:10} \end{equation} Osservato come la seconda equazione del sistema coincida con la condizione \eqref{eq:08}, lo sviluppo delle sommatorie \eqref{eq:10} ci riporta alle equazioni \begin{equation} \cases{\displaystyle \sum x_iy_i= a\sum (x_i)^2+b\sum x_i=0\cr\\ \displaystyle \sum y_i= a \sum x_i+n\,b \cr}\tag{11}\label{eq:eq11} \end{equation} dette "equazioni normali". Queste, dopo alcuni calcoli (aprire, per i particolari, il riquadro sottostante) forniscono i valori cercati per i parametri \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \hat a={n\sum x_iy_i-\bigl(\sum x_i\bigr)\bigl(\sum y_i\bigr)\over n\sum (x_i)^2-\bigl(\sum x_i\bigr)^2},\hskip2cm \hat b={\bigl(\sum y_i\bigr)\sum(x_i)^2-\bigl(\sum x_i y_i\bigr)\sum x_i\over n\sum (x_i)^2-\bigl(\sum x_i\bigr)^2}.\tag{12}\label{eq:12} } \end{equation} La conoscenza di questi due valori permette in definitiva di individuare univocamente l'equazione \(y=\hat a x+\hat b\) rappresentativa della retta che interpola nel modo ottimale le osservazioni sperimentali: possiamo pertanto rappresentarla assieme ai dati così da verificarne, anche visualmente, l'adeguatezza (fig. 4). Tale retta verrà da noi indicata come la retta dei minimi quadrati o, più brevemente, retta di regressione o, ancora, retta di best-fit.
Indicati come valori del best-fit le quantità note \begin{equation} \hat{y}_i=\hat a x_i+\hat b,\qquad i=1,2,\dots, n \end{equation} le differenze \begin{equation} \eqalign{e_i&=y_i-\hat{y}_i\cr &=y_i-\hat a x_i-\hat b,\qquad i=1,2,\dots, n\cr}\tag{13}\label{eq:13} \end{equation} che in figura 4 interpretiamo come le lunghezze con segno dei segmenti verticali, costituiscono i residui. Come già accennato precedentemente (ma nella \eqref{eq:06} rappresentano grandezze ancora incognite), i residui esprimono la discrepanza tra l'osservazione sperimentale e la previsione del modello e, analogamente agli scarti dalla media, la loro somma è nulla come dimostriamo nel riquadro seguente.

Fig. 4. Retta di regressione, dati e rappresentazione dei residui.

Per una verifica della correttezza del modello impiegato che vada oltre un'analisi visiva del grafico di figura 4, risulta utile associare a questo pure il grafico a dispersione dei residui in modo da evidenziare, o meno, loro andamenti anomali (fig. 5).

Fig. 5. Residui della retta \(y=5.19+1.11 x\) di figura 4.
Soluzione delle equazioni normali

Dalla seconda delle \eqref{eq:eq11}, isoliamo a primo membro il parametro \(b\) \begin{equation} b={1\over n}\biggl(\sum y_i-a\sum x_i\biggr)\tag{14}\label{eq:14} \end{equation} e lo sostituiamo nella prima \begin{equation} \sum x_i y_i= a\sum(x_i)^2+{\sum x_i\over n}\biggl(\sum y_i-a\sum x_i\biggr) \end{equation} dalla quale, eseguendo il prodotto, discende \begin{equation} \sum x_i y_i= a\sum(x_i)^2+{\bigl(\sum x_i\bigr)\bigl(\sum y_i\bigr)\over n}-{a\bigl(\sum x_i\bigr)^2\over n}. \end{equation} Possiamo ora risolvere l'equazione nell'unica incognita \(a\) raccogliendola a fattore in un membro \begin{equation} \sum x_i y_i -{\bigl(\sum x_i\bigr)\bigl(\sum y_i\bigr)\over n} = {a\over n}\biggl[n\sum(x_i)^2-\Bigl(\sum x_i\Bigr)^{\kern-3pt 2}\,\biggr] \end{equation} e quindi ottenere la stima di \(a\) \begin{equation} {n\sum x_i y_i -\bigl(\sum x_i\bigr)\bigl(\sum y_i\bigr)\over n} = {a\over n}\biggl[n\sum(x_i)^2-\Bigl(\sum x_i\Bigr)^{\kern-3pt 2}\,\biggr]\qquad\Longrightarrow\qquad \bbox[border:1px solid red,15px,#ffffcc]{\hat{a}={n\sum x_i y_i -\bigl(\sum x_i\bigr)\bigl(\sum y_i\bigr)\over n\sum(x_i)^2-\bigl(\sum x_i\bigr)^2}.}\tag{15}\label{eq:15} \end{equation} Inserito questo risultato nella \eqref{eq:14} riscritta come \begin{equation} n\, b=\sum y_i-\sum x_i\Biggl[{n\sum x_i y_i -\bigl(\sum x_i\bigr)\bigl(\sum y_i\bigr)\over n\sum(x_i)^2-\bigl(\sum x_i\bigr)^2}\Biggr] \end{equation} riportiamo il secondo membro allo stesso denominatore \begin{equation} n\, b={{\sum y_i\biggl[n\sum(x_i)^2-\bigl(\sum x_i\bigr)^2} \biggr]-\sum x_i\biggl[n\sum x_i y_i -\bigl(\sum x_i\bigr)\bigl(\sum y_i\bigr)\biggr]\over n\sum(x_i)^2-\bigl(\sum x_i\bigr)^2} \end{equation} e, a numeratore, svolgiamo i prodotti \begin{equation} n\, b={n\bigl(\sum y_i\bigr)\sum(x_i)^2-\bigl(\sum y_i\bigr)\bigl(\sum x_i\bigr)^2 -n\bigl(\sum x_i\bigr)\bigl(\sum x_i y_i\bigr) +\bigl(\sum x_i\bigr)^2\bigl(\sum y_i\bigr)\over n\sum(x_i)^2-\bigl(\sum x_i\bigr)^2}. \end{equation} Per la presenza di termini opposti, semplifichiamo il numeratore in \begin{equation} n\, b={n\bigl(\sum y_i\bigr)\sum(x_i)^2 -n\bigl(\sum x_i\bigr)\bigl(\sum x_i y_i\bigr) \over n\sum(x_i)^2-\bigl(\sum x_i\bigr)^2} \end{equation} e quindi, dividendo per \(n\), giungiamo al risultato \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \hat b={\bigl(\sum y_i\bigr)\sum(x_i)^2 -\bigl(\sum x_i y_i\bigr)\sum x_i \over n\sum(x_i)^2-\bigl(\sum x_i\bigr)^2}}\tag{16}\label{eq:16} \end{equation}

Somma nulla dei residui

Ripresa la definizione di residuo \eqref{eq:13} \begin{equation} e_i=y_i-\hat a x_i-\hat b,\qquad i=1,2,\dots,n \end{equation} con \(\hat a\) e \(\hat b\) espresse dalle \eqref{eq:12}, consideriamo la loro somma \begin{equation} \sum e_i= \sum(y_i-\hat a x_i-\hat b)\qquad\Longrightarrow\qquad \sum e_i=\sum y_i-\hat a\sum x_i-n\, \hat b. \end{equation} Posto per brevità \begin{equation} \Delta= n\sum (x_i)^2-\Bigl(\sum x_i\Bigr)^2 \end{equation} inseriamo nella precedente i risultati \eqref{eq:12} \begin{equation} \sum e_i= \sum y_i-\sum x_i\Biggl[{n\sum x_iy_i-\bigl(\sum x_i\bigr)\bigl(\sum y_i\bigr)\over \Delta}\Biggr]-n \Biggl[{\bigl(\sum y_i\bigr)\sum(x_i)^2-\bigl(\sum x_i y_i\bigr)\sum x_i\over \Delta}\Biggr] \end{equation} e riportiamo il secondo membro ad un comun denominatore \begin{equation} \sum e_i= {1\over \Delta}\sum y_i\Bigl[n\sum (x_i)^2-\Bigl(\sum x_i\Bigr)^2\Bigr]-\sum x_i\Bigl[n\sum x_iy_i-\Bigl(\sum x_i\Bigr)\Bigl(\sum y_i\Bigr)\Bigr]-n \Bigl[\Bigl(\sum y_i\Bigr)\sum(x_i)^2-\Bigl(\sum x_i y_i\Bigr)\sum x_i\Bigr]. \end{equation} Svolti i prodotti tra le sommatorie \begin{equation} \sum e_i= {1\over \Delta}\Bigl[n\Bigl(\sum y_i\Bigr)\sum (x_i)^2-\Bigl(\sum y_i\Bigr)\Bigl(\sum x_i\Bigr)^2-n\sum x_i\Bigl(\sum x_iy_i\Bigr)+\Bigl(\sum x_i\Bigr)^2\Bigl(\sum y_i\Bigr)-n \Bigl(\sum y_i\Bigr)\sum(x_i)^2+n\Bigl(\sum x_i y_i\Bigr)\sum x_i\Bigr]. \end{equation} notiamo la presenza di tre coppie di termini opposti, per cui rimane dimostrato che \begin{equation} \sum e_i= {1\over \Delta}\bigl[\,0\,\bigr]=0. \end{equation}

Caso particolare: la proporzionalità

Un caso particolare di regressione lineare si presenta quando il modello che lega le grandezze \(X\) e \(Y\) si riduce a quello di una semplice proporzionalità ossia \begin{equation} y=a x.\tag{17}\label{eq:17} \end{equation} In tal caso la retta che lo rappresenta possiede intercetta nulla e passa, obbligatoriamente, per l'origine. In questo modello la funzione da minimizzare \eqref{eq:09} diviene \begin{equation} s^2(a)=\sum_{i\,=\,1}^n e_i^2=\sum_{i\,=\,1}^n (y_i-a x_i)^2.\tag{18}\label{eq:18} \end{equation} e quindi appare dipendente da una sola variabile. La sua derivata è quindi \begin{equation} {d(s^2)\over dt}=\sum_{i\,=\,1}^n \Bigl[-2x_i(y_i-a x_i)\Bigr] \end{equation} per cui, posta questa a zero, \begin{equation} {d(s^2)\over dt}=0\qquad\Longrightarrow\qquad \sum_{i\,=\,1}^n \Bigl[-2x_i(y_i-a x_i)\Bigr]=0 \end{equation} il valore \(\hat a\) di best fit cui corrisponde il minimo di \(s^2\) è \begin{equation} \sum_{i\,=\,1}^n (-x_iy_i)+a \sum_{i\,=\,1}^n(x_i)^2=0 \end{equation} \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \hat{a}={\sum x_iy_i\over \sum (x_i)^2}\tag{19}\label{eq:19} } \end{equation} Ottenuto questo valore, i residui assumono la forma \begin{equation} e_i=y_i-\hat{y}_i=y_i-\hat{a}x_i\qquad i=1,2,\dots, n.\tag{20}\label{eq:20} \end{equation}

A discapito comunque della semplicità della relazione \eqref{eq:19}, l'imposizione del passaggio per l'origine della retta che rappresenta il modello, implica non trascurabili condizionamenti sui residui che, come mostreremo nella pagina degli esempi di laboratorio (v. legge di Ohm), non possiedono necessariamente somma nulla. Ne segue che la loro deviazione standard \(\sigma\) è, in genere, maggiore della corrispondente deviazione nel modello lineare completo. Inoltre alcuni importanti indici, quali il coefficiente di correlazione discusso in una successiva sezione e rappresentato dalla relazione \eqref{eq:35} perdono la loro validità e vanno sostituiti con altre espressioni.

Ipotesi sulle incertezze dei parametri

Nella sezione precedente, definita una funzione lineare, abbiamo affrontato il problema della ricerca di quella particolare relazione in grado di descrivere al meglio il legame tra due grandezze \(X\) e \(Y\) per le quali si dispone di una serie di coppie di dati sperimentali. Definito il criterio dei minimi quadrati, è stato poi possibile determinare la coppia di valori, \(\hat a\) e \(\hat b\), che individuano la pendenza e l'intercetta della retta ottimale. Come accennato inizialmente non abbiamo invece considerata la presenza di errori e incertezze che invece pure caratterizzano ciascuna misura. In questa sezione intendiamo perciò studiare come queste incertezze si riflettano sui valori dei parametri rappresentati dalle \eqref{eq:12}.

Le ipotesi che poniamo per poter trattare almeno il caso più frequente sono:

  • la variabile \(x\) è priva di errore oppure, data la piccolezza di questo, è possibile, con buona approssimazione, trascurare la sua incertezza.
  • la medesima incertezza \(\sigma\) caratterizza tutti i valori \(y_i\) della variabile \(y\) e questa è assunta normalmente distribuita attorno al valore vero.

Se quindi, per es. in base a considerazioni collegate agli strumenti utilizzati, supponiamo noto l'errore \(\sigma\) della variabile dipendente, diviene possibile rendere il grafico a dispersione di figura 4 più realistico associando ad ogni misura \(y_i\) il rispettivo intervallo di incertezza \([y_i-\sigma, y_i+\sigma]\) (fig. 6): in tal modo appare con maggiore evidenza la necessità di definire le incertezze sui parametri della retta di regressione.

Fig. 6. Retta di regressione, dati e ipotetici errori.

Propagazione degli errori e stima delle incertezze dei parametri

Per determinare l'incertezza su \(\hat a\) e \(\hat b\), dobbiamo innanzitutto considerare come si calcola la deviazione standard \(\sigma_z\) di una grandezza \(Z\) dipendente da un certo numero di parametri \(y_i\) con \(i=1\dots n\), ciascuno dei quali possiede la sua propria deviazione standard \(\sigma_i\). La formula che ci permette questo calcolo è \begin{equation} \sigma_z^2=\sum_{i\,=\,1}^n\Bigl({\partial z\over \partial y_i}\sigma_i\Bigr)^2, \tag{21}\label{eq:21} \end{equation} cosicché la deviazione standard cercata \(\sigma_z\) è la radice quadrata del valore dato dall'espressione \eqref{eq:21}.
Se quindi consideriamo le espressioni \eqref{eq:12} per \(\hat a\) e \(\hat b\) come funzioni delle \(y_i\), i quadrati delle rispettive deviazioni standard sono \begin{equation} \sigma_a^2=\sum_{i\,=\,1}^n\Bigl({\partial a\over \partial y_i}\sigma_i\Bigr)^2,\hskip2cm \sigma_b^2=\sum_{i\,=\,1}^n\Bigl({\partial b\over \partial y_i}\sigma_i\Bigr)^2 \end{equation} e sfruttando l'ipotesi che suppone tutte uguali la singole deviazioni cioè \(\forall\, i=1\dots n,\ \sigma_i=\sigma\), queste si semplificano nelle \begin{equation} \sigma_a^2=\sigma^2\sum_{i\,=\,1}^n\Bigl({\partial a\over \partial y_i}\Bigr)^2\hskip2cm \sigma_b^2=\sigma^2\sum_{i\,=\,1}^n\Bigl({\partial b\over \partial y_i}\Bigr)^2. \tag{22}\label{eq:22} \end{equation} Il calcolo esplicito riportato nel riquadro seguente fornisce i risultati \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \sigma_a=\sigma\sqrt{n\over \Delta},\hskip2cm \sigma_b=\sigma\sqrt{\sum (x_i)^2\over \Delta}\qquad\hbox{con}\qquad \Delta=n\sum(x_i)^2-\Bigl(\sum x_i\Bigr)^2 }\tag{23}\label{eq:23} \end{equation} evidentemente dipendenti dalla comune deviazione standard \(\sigma\). Per poter quindi utilizzare queste ultime relazioni è necessario fornire comunque una stima della deviazione \(\sigma\) delle misure \(y_i\). Poiché queste sono state ottenute sperimentalmente e, presumibilmente, per mezzo di un qualche strumento, potremmo assegnare a \(\sigma\) l'errore massimo che caratterizza lo strumento. Se invece non disponiamo di una tale informazione, non ci rimane che calcolare la deviazione \(\sigma\) a partire da una "media" dei quadrati dei residui. Ripresa quindi la definizione di residuo \eqref{eq:13} la nostra stima della deviazione standard complessiva è data dalla formula \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \sigma^2={1\over n-2}\sum_{i\,=\,1}^n e_i^2 ={1\over n-2}\sum_{i\,=\,1}^n\bigl(y_i-\hat a x_i-\hat b\bigr)^2 } \tag{24}\label{eq:24} \end{equation} dove il fattore \(n-2\) sta in luogo di \(n\) in quanto rappresenta l'effettivo numero di grandezze indipendenti essendo \(\hat a\) e \(\hat b\) valori noti. Questa quantità è anche detta media dei quadrati dei residui mentre \(\sigma\) è chiamato errore standard della regressione e poiché questa grandezza è calcolata a partire dai residui che, evidentemente dipendono dalla funzione di fit \eqref{eq:05}, il suo valore è sensibilmente influenzato dal modello scelto piuttosto che da fattori esterni.
Con queste ultime osservazioni siamo in grado di associare alla retta di regressione nella figura 7 anche la sua fascia di incertezza: questa è delimitata inferiormente e superiormente dalle due rette \[ y_{inf}=(a-\sigma_a)x+(b-\sigma_b),\qquad y_{sup}=(a+\sigma_a)x+(b+\sigma_b). \] Nella figura 8 aggiungiamo alla già riportata distribuzione dei residui pure la fascia che comprende i valori interni all'intervallo \([-\sigma,+\sigma]\) con \(\sigma\) errore standard della regressione.

Fig. 7. Retta di regressione e sua fascia di incertezza.
Fig. 8. Residui e deviazione standard \(\sigma\).
Dimostrazione delle deviazioni standard di \(a\) e \(b\)

La dimostrazione delle formule \eqref{eq:23} parte innanzitutto con il calcolo delle derivate parziali ad argomento delle sommatorie \eqref{eq:22}. Considerando che ciascuna variabile \(y_i\) nella sommatoria appare a primo grado, il calcolo delle due derivate parziali fornisce \begin{equation} {\partial a\over \partial y_i}={n x_i-\sum x_i\over \Delta},\hskip2cm {\partial b\over \partial y_i}={\sum(x_i)^2-x_i\sum x_i\over \Delta}.\tag{25}\label{eq:25} \end{equation} Svolgendo il quadrato della derivata di \(a\) \begin{equation} \eqalign{\Bigl({\partial a\over \partial y_i}\Bigr)^2&=\Bigl({n x_i-\sum x_i\over \Delta}\Bigr)^2\cr &={n^2x_i^2+(\sum x_i)^2-2n x_i\sum x_i\over \Delta^2},\cr} \end{equation} e sommando sull'indice \(i\) di questa espressione, otteniamo come risultato \begin{equation} \eqalign{\sum_{i\,=\,1}^n\Bigl({\partial a\over \partial y_i}\Bigr)^2&={1\over \Delta^2}\biggl[{n^2\sum (x_i)^2+n\Bigl(\sum x_i\Bigr)^2-2n\Bigl(\sum x_i\Bigr)^2}\biggr]\cr &={n\over\Delta^2}\biggl[n \sum(x_i)^2-\Bigl(\sum x_i\Bigr)^2\biggr]={n\Delta\over \Delta^2}={n\over\Delta}.\cr} \end{equation} Associato questo termine ai quadrati delle deviazioni, \begin{equation} \sigma_a^2=\sigma^2\cdot {n\over \Delta} \end{equation} l'estrazione della radice quadrata ci fornisce la deviazione standard cercata \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \sigma_a=\sigma\sqrt{n\over\Delta}. } \end{equation} Procediamo allo stesso modo per \(\partial b/\partial y_i\). Il quadrato della seconda delle \eqref{eq:25} è l'espressione \begin{equation} \eqalign{\Bigl({\partial b\over \partial y_i}\Bigr)^2&=\Biggl[{\sum(x_i)^2-x_i\sum x_i\over \Delta}\Biggr]^2\cr &={\Bigl[\sum(x_i)^2\Bigr]^2+x_i^2\Bigl(\sum x_i\Bigr)^2-2x_i\Bigl(\sum x_i\Bigr)\sum(x_i)^2\over \Delta^2},\cr} \end{equation} che inserita nella sommatoria, riscriviamo quest'ultima come \begin{equation} \eqalign{\sum_{i\,=\,1}^n\Bigl({\partial b\over \partial y_i}\Bigr)^2 &={1\over\Delta^2} \biggl\{ n\Bigl[\sum(x_i)^2\Bigr]^2 +\sum (x_i)^2\Bigl(\sum x_i\Bigr)^2 -2\sum(x_i)^2\Bigl(\sum x_i\Bigr)^2 \biggr\}\cr &={1\over\Delta^2}\biggl\{n\Bigl[\sum(x_i)^2\Bigr]^2-\sum(x_i)^2\Bigl(\sum x_i\Bigr)^2\biggr\}.\cr} \end{equation} Posto a fattore il termine comune, la sommatoria si semplifica in \begin{equation} \eqalign{\sum_{i\,=\,1}^n\Bigl({\partial b\over \partial y_i}\Bigr)^2 &={1\over\Delta^2}\Bigl\{{\sum(x_i)^2\Bigl[{n\sum(x_i)^2-\Bigl(\sum x_i\Bigr)^2}\Bigr]}\Bigr\} \cr &={1\over\Delta^2}\Bigl\{\sum (x_i)^2\cdot \Delta\Bigr\}={1\over\Delta}{\sum (x_i)^2}.\cr} \end{equation} Il quadrato della deviazione \(\sigma_b\) è pertanto \begin{equation} \sigma_b^2=\sigma^2\cdot {1\over\Delta}{\sum (x_i)^2}, \end{equation} dalla quale, con l'estrazione della radice quadrata, troviamo \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \sigma_b=\sigma\sqrt{\sum(x_i)^2\over \Delta}. } \end{equation}

Caso particolare: incertezza nel parametro della proporzionalità

Il calcolo dell'incertezza del parametro \(\hat{a}\) nel modello con intercetta nulla o della proporzionalità segue le stesse orme di quello generale a partire dalla relazione \eqref{eq:21} \begin{equation} \sigma_a^2=\sum_{i\,=\,1}^n\Bigl({da\over d y_i}\sigma_i\Bigr)^2 \tag{26}\label{eq:26} \end{equation} che, nell'ipotesi di uguali deviazioni standard \(\sigma_i=\sigma,\quad i=1,\,\dots,\,n\) per le misure \(y_i\) si riduce alla \begin{equation} \sigma_a^2=\sigma^2\sum_{i\,=\,1}^n\Bigl({d a\over d y_i}\Bigr)^2. \tag{27}\label{eq:27} \end{equation} Poiché, in base alla \eqref{eq:19} \begin{equation} {da\over dy_i}={x_i\over \sum(x_i)^2}\quad\Longrightarrow\quad \Bigl({da\over dy_i}\Bigr)^2={x_i^2\over \bigl[\sum(x_i)^2\bigr]^2}, \end{equation} la \eqref{eq:27} diviene \begin{equation} \sigma_a^2=\sigma^2\sum_{i\,=\,1}^n\Bigl({da\over dy_i}\Bigr)^2={\sigma^2 \sum (x_i)^2\over \bigl[\sum(x_i)^2\bigr]^2}={\sigma^2\over \sum(x_i)^2}\tag{28}\label{eq:28} \end{equation} dalla quale otteniamo \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \sigma_a={\sigma\over \sqrt{\sum(x_i)^2}}\tag{29}\label{eq:29} } \end{equation} dove la deviazione \(\sigma\) dei residui, \eqref{eq:24}, è espressa dalla \begin{equation} \sigma^2={1\over n-1}\sum_{i\,=\,1}^n e_i^2 ={1\over n-1}\sum_{i\,=\,1}^n\bigl(y_i-\hat a x_i\bigr)^2 \end{equation} e quindi \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ \sigma=\sqrt{{\sum (y_i-\hat a x_i)^2\over n-1}} } \tag{30}\label{eq:30} \end{equation} in quanto dovendo, in questo caso, determinare il solo parametro \(a\), sono \(n-1\) i gradi di libertà.

Coefficiente di correlazione

L'insieme \(A\) dei dati iniziali consiste di coppie di misure \((x_i,y_i)\) dove la quantità \(y\) è supposta dipendere dalla variabile \(x\). L'ipotesi finora seguita consiste nell'affermare l'esistenza di una relazione lineare della forma \(y=ax+b\) tra le variabili \(x\) e \(y\) e su tali basi abbiamo sviluppato nelle sezioni precedenti un percorso che ci ha permesso di determinare l'intercetta \(\hat b\) e la pendenza \(\hat a\) della retta che meglio interpola questi dati sperimentali. In questa sezione intendiamo definire una qualche grandezza che ci possa informare in quale misura la retta ottimale si adatti ai nostri dati di partenza in tal modo fornendo una conferma o al contrario, una smentita, al modello lineare scelto.

Sappiamo che due grandezze sono indipendenti se al variare dell'una, per es. \(x\), l'altra rimane costante ossia, poste entrambe in un grafico, i valori della seconda dovrebbero distribuirsi su una retta orizzontale. Sulla base di tale osservazione potremmo quindi analizzare il valore del coefficiente \(\hat a\) di \eqref{eq:12} e, nell'ipotesi di indipendenza, aspettarci che esso sia nullo e quindi che la retta ottimale sia orizzontale e parallela all'asse della variabile \(x\). D'altra parte, a causa della dispersione dei dati e quindi dell'errore \(\sigma_a\), potremmo trovare una pendenza prossima allo zero e, per questo, non essere ancora in grado di discriminare tra la dipendenza o l'indipendenza delle quantità coinvolte. Poiché la dispersione dei dati è sostanzialmente sempre presente, il parametro \(\hat a\) non sembra pertanto adatto a tale scopo per cui appare necessario ricercare un'alternativa in grado di fornirci informazioni più precise sulla relazione che lega le due grandezze.

A tal fine calcoliamo le medie dei valori di \(x\) e \(y\) \begin{equation} \overline{x}={1\over n}{\sum_{i\,=\,1}^n x_i},\hskip2cm \overline{y}={1\over n}{\sum_{i\,=\,1}^ny_i}\tag{31}\label{eq:31} \end{equation} e nel grafico a dispersione dei dati aggiungiamo la retta verticale \(x=\overline{x}\) e la retta orizzontale \(y=\overline{y}\) (fig. 9)

Fig. 9. Retta di regressione e segno dei prodotti \((x_i-\overline{x})(y-\overline{y}))\).

Queste linee dividono il grafico in quattro quadranti e per ogni coppia \((x_i,y_i)\) calcoliamo le quantità

  • \(x_i-\overline{x}\), cioè la deviazione della misura \(x_i\) dalla media della variabile \(x\),
  • \(y_i-\overline{y}\), cioè la deviazione della misura \(y_i\) dalla media della variabile \(y\),
  • il prodotto delle precedenti quantità \((x_i-\overline{x})(y_i-\overline{y})\).

È evidente che la grandezza \(x_i-\overline{x}\) sarà positiva nel I e IV quadrante e negativa nel II e III, mentre \(y_i-\overline{y} > 0\) nel I e II e \(y_i-\overline{y} < 0\) nel III e IV quadrante (tabella 1).

quadrante
\(x_i-\overline{x}\)
\(y_i-\overline{y}\)
\((x_i-\overline{x})(y_i-\overline{y})\)
1
\(+\)
\(+\)
\(+\)
2
\(-\)
\(+\)
\(-\)
3
\(-\)
\(-\)
\(+\)
4
\(+\)
\(-\)
\(-\)
 
Tabella 1. Segno delle quantità \(x_i-\overline{x}\), \(y_i-\overline{y}\) e del loro prodotto.

Uno sguardo alla figura 9 mostra che, se la relazione tra \(x\) e \(y\) ha pendenza positiva, allora vi è un numero maggiore di punti nel I e III quadrante piuttosto che nel II e IV mentre avviene il contrario se la pendenza fosse negativa (fig. 10)

Fig. 10. Retta di regressione e segno dei prodotti \((x_i-\overline{x})(y-\overline{y}))\).

Nel primo caso la somma dei prodotti riportati nella quarta colonna della tabella 1 sarà molto probabilmente positiva, mentre l'opposto avviene per una relazione con pendenza negativa.
Ne segue che se consideriamo il valor medio della somma dei prodotti \((x_i-\overline{x})(y-\overline{y}))\), il segno di questa nuova grandezza detta covarianza tra le quantità \(X\) e \(Y\) (più avanti chiariremo il perché del fattore \(n-1\)) \begin{equation} \hbox{cov}(X,Y)={1\over n-1}\sum_{i\,=\,1}^n (x_i-\overline{x})(y_i-\overline{y})\tag{32}\label{eq:32} \end{equation} ci informa sulla pendenza della relazione funzionale \(y=f(x)\) esistente tra le variabili \(x\) e \(y\):

  • se \(\hbox{cov}(X,Y)> 0\) allora la monotonia di \(y=f(x)\) è crescente mentre,
  • se \(\hbox{cov}(X,Y)< 0\) allora la monotonia di \(y=f(x)\) è decrescente.

La covarianza di due grandezze è evidentemente una quantità espressa dal prodotto delle unità di misura scelte per \(X\) e per \(Y\). Ne segue che il suo valore numerico dipende da queste unità e può cambiare a seconda delle unità utilizzate. Per disporre di una informazione che sia invece indipendente dalle unità conviene definire delle grandezze che siano adimensionali. Pertanto, ricordando che ad un valor medio si affianca in genere la sua deviazione standard o scarto quadratico medio, possiamo standardizzare le variabili originarie sostituendole nella \eqref{eq:32} con \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ r(X,Y)={1\over n-1}\sum_{i\,=\,1}^n\Bigl({x_i-\overline{x}\over \sigma_x}\Bigr)\Bigl({y_i-\overline{y}\over \sigma_y}\Bigr)\tag{33}\label{eq:33} } \end{equation} dove \begin{equation} \sigma_x=\sqrt{\sum(x_i-\overline{x})^2\over n-1},\hskip2cm \sigma_y=\sqrt{\sum (y_i-\overline{y})^2\over n-1}\tag{34}\label{eq:34} \end{equation} sono rispettivamente le deviazioni standard di \(X\) e \(Y\).

Questa nuova grandezza, detta coefficiente di correlazione tra \(X\) e \(Y\) conserva le caratteristiche del segno della covarianza ma, come richiesto, è indipendente dalle unità. Poiché inoltre nella stima degli scarti quadratici medi \eqref{eq:34} appare a denominatore il termine \(n-1\), la presenza del fattore \(1/(n-1)\) nel coefficiente di correlazione ci permette di riscriverlo come \begin{equation} \bbox[border:1px solid red,15px,#ffffcc]{ r(X,Y)={ \sum (x_i-\overline{x})(y_i-\overline{y})\over \sqrt{\sum (x_i-\overline{x})^2\sum (y_i-\overline{y})^2} }\tag{35}\label{eq:35} } \end{equation} in quanto il prodotto delle deviazioni è \begin{equation} \sigma_x\sigma_y= {1\over n-1}\sqrt{\sum (x_i-\overline{x})^2\sum (y_i-\overline{y})^2}. \end{equation} A seguito di questa standardizzazione delle variabili sostituiamo i precedenti grafici a dispersione 9 e 10, rispettivamente, con quelli di fig. 11 e 12.

Fig. 11. Segno dei prodotti \((x_i-\overline{x})(y-\overline{y}))/(\sigma_x\sigma_y)\).
Fig. 12. Segno dei prodotti \((x_i-\overline{x})(y-\overline{y}))/(\sigma_x\sigma_y)\)\).

Una proprietà importante del coefficiente di correlazione

In aggiunta al suo segno, una fondamentale proprietà del coefficiente di correlazione è rappresentata dal valore numerico che esso assume al variare della dispersione dei dati attorno alla retta di regressione. Difatti osservando le tre coppie di figure sottostanti appare evidente che, al diminuire della dispersione dei dati, il valore assoluto da questo coefficiente tenda ad 1 e lo raggiunga quando gli stessi dati appartengano tutti alla medesima retta ossia abbiano una dispersione nulla.

Fig. 13. Dispersione dei dati e coefficiente di correlazione.
Fig. 14. Dispersione dei dati e coefficiente di correlazione.
Fig. 15. Dispersione nulla dei dati e coefficiente di correlazione unitario.

Quest'ultima notevole proprietà si dimostra facilmente in quanto, nell'ipotesi che tutti i punti \((x_i,y_i)\) appartengano alla retta di equazione \(y=\hat ax+\hat b\) allora, assieme alla \begin{equation} y_i=\hat a x_i+\hat b,\tag{36}\label{eq:36} \end{equation} vale pure la relazione che coinvolge la somma \begin{equation} {\sum y_i}=\hat a{\sum x_i}+\sum b\qquad\Longrightarrow\qquad {\sum y_i}=\hat a{\sum x_i}+n\, \hat b\qquad\Longrightarrow\qquad {1\over n}{\sum y_i}=\hat a\Bigl({1\over n}{\sum x_i}\Bigr)+\hat b \end{equation} dalla quale discende pure il legame tra i valori medi \begin{equation} \overline{y}=\hat a\,\overline{x}+\hat b.\tag{37}\label{eq:37} \end{equation} Sottraendo quest'ultima dalla \eqref{eq:36} \begin{equation} y_i-\overline{y}=a (x_i-\overline{x})\tag{38}\label{eq:38} \end{equation} e inserita nella espressione \eqref{eq:35} del coefficiente di correlazione, questo risulta \begin{equation} r(X,Y)={\sum(x_i-\overline{x})a(x_i-\overline{x})\over \sqrt{\sum(x_i-\overline{x})^2a^2\sum(x_i-\overline{x})^2}}.\tag{39}\label{eq:39} \end{equation} Eseguendo i prodotti e semplificando, otteniamo la tesi \begin{equation} \eqalign{r(X,Y)&= {a\sum (x_i-\overline{x})^2\over \sqrt{a^2[\sum(x_i-\overline{x})^2]^2}}\cr &={a\sum (x_i-\overline{x})^2\over |a| \sum (x_i-\overline{x})^2}\cr &={a\over |a|}=\pm 1.\cr }\tag{40}\label{eq:40} \end{equation} In generale comunque, poiché i punti sperimentali sono caratterizzati da una, seppur piccola, dispersione, il coefficiente di correlazione soddisfa alle disuguaglianze \begin{equation} -1≤r(X,Y)≤+1. \end{equation} Infine, nel caso che i dati si distribuiscano casualmente senza suggerire una qualche relazione tra \(X\) e \(Y\), ci aspettiamo che il valore del coefficiente di correlazione sia prossimo allo zero come, in effetti, mostra la figura 16.

Fig. 16. Assenza di correlazione tra le grandezze \(X\) e \(Y\).

In definitiva, il coefficiente di correlazione \(r(X,Y)\) misura il grado di correlazione lineare tra due variabili: esso è nullo quando non c'è alcuna correlazione (lineare!) mentre vale \(\pm 1\) se la correlazione è completa.

Modelli alternativi e adattamento ad altre curve

Il modello lineare seguito finora si basa su una relazione del tipo \(y=ax+b\) tra le variabili \(x\) e \(y\) rappresentative delle grandezze \(X\) e \(Y\). Ovviamente, come possiamo rilevare dalla figura 17 dove, ad un nuovo insieme di dati, associamo il grafico della retta di regressione e a lato l'andamento dei residui con la loro deviazione standard \(\sigma\), questo modello non può rappresentare l'unica scelta.

Fig. 17. Necessità di un modello alternativo tra le grandezze \(X\) e \(Y\).

Già a livello visuale, la distribuzione delle misure in figura 17, evidentemente lontana dal suggerire il modello lineare, rende la ricerca di un legame alternativo tra \(X\) e \(Y\) particolarmente difficile in quanto molte sono le possibili curve che si possono adattare ai rilievi sperimentali. Per tale motivo e come accennato inizialmente, è spesso conveniente disporre di qualche informazione aggiuntiva di carattere teorico o sperimentale che sia di supporto alla scelta del modello. In alcuni casi comunque, la ricerca di un modello non lineare si può ricondurre ai metodi sviluppati per quello lineare e, in questa sezione, presentiamo due semplici esempi.

Modello esponenziale. Supponiamo quindi che le misure di figura 17 siano descritte da una funzione esponenziale del tipo \begin{equation} y=a\,e^{b\,x},\tag{41}\label{eq:41} \end{equation} funzione peraltro molto frequente in ambito fisico. Questa, come nel modello lineare, dipende dai parametri \(a\) e \(b\) e di questi intendiamo ottenere delle stime in modo che il grafico di \eqref{eq:41} si adatti nel modo migliore alla distribuzione dei dati (fig. 17). A partire dal modello rappresentato dalla \eqref{eq:41}, possiamo riscrivere questa funzione considerando il logaritmo naturale di entrambi i membri: otteniamo \begin{equation} \log{y}=\log\bigl(a e^{bx}\bigr)\quad\Longrightarrow\quad \log{y}=\log{a}+b x\tag{42}\label{eq:42} \end{equation} per cui, introdotta la nuova variabile \(z=\log{y}\), la relazione tra quest'ultima e la variabile \(x\) assume la forma lineare \begin{equation} z=\log{a}+b\, x.\tag{43}\label{eq:43} \end{equation} Ne segue che se disponiamo sull'asse verticale di un nuovo grafico a dispersione i valori \(z_i=\log{y_i}\), la distribuzione delle coppie \((x_i,z_i)\) dovrebbe suggerire la correttezza dell'ipotesi \eqref{eq:41} oppure, eventualmente, smentirla. Nel caso in esame i dati così trasformati (fig. 18)

Fig. 18. Modello lineare tra le grandezze \(X\) e \(Z=\log(Y)\).

confermano in modo abbastanza evidente che il legame ipotizzato sembra essere corretto. Si giustifica pertanto il calcolo della retta di regressione ottimale \(z=\hat c+\hat b\, x\) dove i parametri \(\hat b\) e \(\hat c\) sono calcolati applicando le formule \eqref{eq:12}. Dalla conoscenza di \(\hat c\) e tramite la trasformazione inversa del logaritmo, otteniamo la stima del parametro originario \(\hat a\) \[ \hat c= \log{\hat a}\qquad\Longrightarrow\qquad \hat a= \exp({\hat c}). \] In tal modo la funzione esponenziale ottimale cercata è \[ y=\hat a\, e^{\hat b\, x}. \] Nella figura 19 sovrapponiamo ai dati di partenza la curva esponenziale ottimale e, a fianco, riportiamo la distribuzione dei residui con in evidenza l'intervallo \([-\sigma,+\sigma]\) dove la deviazione standard \(\sigma\) è calcolata in base alla formula \begin{equation} \sigma^2={1\over n-2}\sum_{i\,=\,1}^n e_i^2 ={1\over n-2}\sum_{i\,=\,1}^n\Bigl(y_i-\hat a\, e^{\hat b\, x_i}\Bigr)^2 \tag{44}\label{eq:44} \end{equation} ovvia generalizzazione della \eqref{eq:24}. L'assenza di un qualche andamento sistematico dei residui è un ulteriore indice della correttezza dell'ipotesi \eqref{eq:41} su cui si basa il modello. Le incertezze \(\sigma_{\hat c}\) e \(\sigma_{\hat b}\) sui parametri discendono invece dalla applicazione delle relazioni \eqref{eq:23}: infine \(\sigma_{\hat a}\) deriva dalla \(\sigma_{\hat c}\) prendendone l'esponenziale.

Fig. 19. Curva esponenziale di regressione tra \(X\) e \(Y\) e residui.

Modello omografico. In questo secondo esempio, disponiamo dei dati rappresentati nel diagramma di figura 20 e con questa base sperimentale riportiamo la retta di regressione la cui pendenza e intercetta si sono ottenute utilizzando le formule \eqref{eq:12}. Evidentemente la retta e l'andamento sistematico dei residui mostrano con chiarezza l'inadeguatezza del modello lineare.

Fig. 20. Necessità di un modello non lineare tra \(X\) e \(Y\).

Facciamo quindi l'ipotesi che il modello seguito dalle grandezze \(X\) e \(Y\) sia invece descritto dalla funzione non lineare \begin{equation} y= {a\over x}+b\tag{45}\label{eq:45} \end{equation} dipendente dai due parametri \(a\) e \(b\). Questa funzione rappresenta un caso particolare della funzione omografica (da qui il nome convenzionale che abbiamo dato al modello) \begin{equation} y= {a+bx\over cx+d}\tag{46}\label{eq:46} \end{equation} il cui grafico è un'iperbole equilatera traslata. La funzione \eqref{eq:45}, se riscritta come \((y-b)x=a\), può inoltre rappresentare la proporzionalità inversa tra le grandezze \(X\) e \(Y-b\). Comunque, nell'ipotesi che sia \(a>0\), questo modello assicura una monotonia decrescente della variabile \(y\) al crescere della variabile positiva \(x\) e quindi appare coerente con l'andamento delle misure.
Diversamente dal caso precedente dove, per ottenere la linearizzazione abbiamo definito la nuova grandezza \(Z\) trasformando la \(Y\), qui definiamo la nuova \(Z\) a partire dalla \(X\) ponendo tra le rispettive variabili la relazione \begin{equation} z={1\over x}.\tag{47}\label{eq:47} \end{equation} Ciò permette di riscrivere la \eqref{eq:45} nella forma lineare \begin{equation} y=a z+b\tag{48}\label{eq:48} \end{equation} dove i parametri \(a\) e \(b\) sono gli stessi della forma non lineare. In base a ciò, la distribuzione delle coppie \((z_i,y_i)\) con \(z_i=1/x_i\) è ora quella riportata in figura 21 e pur con gli usuali scarti dall'andamento lineare, appare giustificato il calcolo della retta di regressione \(y=\hat a\,z+\hat b\) tra le grandezze \(Z=1/X\) e \(Y\).

Fig. 21. Modello lineare tra le grandezze \(Z=1/X\) e \(Y\).

Per verificarne la compatibilità con le coppie iniziali \((x_i,y_i)\) non è necessaria alcune trasformazione dei parametri \(\hat a\) e \(\hat b\) e quindi possiamo sovrapporre il grafico della funzione ottimale \begin{equation} y={\hat a\over x}+\hat b\tag{49}\label{eq:49} \end{equation} a quello originario (fig. 22). La distribuzione casuale dei residui conferma ulteriormente la bontà della scelta iniziale del modello omografico. Naturalmente le incertezze discendono dalla \eqref{eq:23} con l'accorgimento di sostituire le \(x_i\) con \(z_i=1/x_i\).

Fig. 22. Dipendenza omografica tra \(X\) e \(Y\) e residui.

Concludiamo questa sezione con un confronto tra modelli alternativi ma a partire dalla stessa base di dati.
Facciamo quindi l'ipotesi che il modello alternativo a quello omografico sia quello esponenziale e procediamo, come mostrato nel primo esempio, al calcolo dei parametri della funzione \[ y=\hat a\, e^{\hat b\,x} \] e nelle figure 23 e 24 riportiamo graficamente gli esiti delle elaborazioni.

Fig. 23. Modello esponenziale tra \(X\) e \(Y\) e residui.

La linearizzazione delle coppie di misure originarie \((x_i,y_i)\) con \(i=1,2,\dots, n\) ottenuta con la trasformazione delle ordinate \(z_i=\log(y_i)\) avvicina l'insieme ad una retta ma i risultati non appaiono del tutto soddisfacenti. In entrambi i grafici di figura 23 le curve "ottimali" tracciate non si adattano ai dati e questi non si distribuiscono in modo casuale attorno a queste. Tutto ciò risulta più evidente nel grafico dei residui (fig. 24) dove l'andamento complessivo possiede un evidente carattere sistematico che comporta, nel caso del modello esponenziale, una deviazione standard pari a \(\sigma=0.39\) sensibilmente maggiore di quella del modello omografico \(\sigma=0.14\).

Fig. 24. Residui nel modello esponenziale.

In definitiva, tra i due modelli, quello omografico appare più aderente ai dati sperimentali.

Per altri esempi di linearizzazione si veda la pagina di laboratorio in corrispondenza delle esperienze su: Moto uniformemente accelerato, Legge di Boyle, Legge di Snell.

PARTE 2

Generalizzazione del metodo dei minimi quadrati

Il percorso seguito finora nella ricerca di un modello in grado di sintetizzare un insieme di osservazioni sperimentali si è focalizzato essenzialmente sulla sola funzione lineare \(f(x)=a x+b\) o, al più, su funzioni riconducili a questa. Qui, sulla base del diagramma a blocchi di figura 25, intendiamo mostrare come il metodo dei minimi quadrati si possa generalizzare così da comprendere un più ampio spettro di funzioni.

Fig. 25. Sviluppo del metodo dei minimi quadrati.
  1. Come più volte evidenziato, sulla base di solide considerazioni fisiche di carattere teorico o sperimentali, definiamo innanzitutto la funzione modello o funzione di fit \begin{equation} y=f(x,a,b,\dots).\tag{50}\label{eq:50} \end{equation} Questa funzione dipende, nei casi più semplici, da una variabile indipendente \(x\) e da un certo numero di parametri incogniti \(a,\,b,\,\dots\).
  2. Avendo a disposizione l'insieme di coppie \((x_i,y_i)\) con \(i=1,\,2\,\dots,n\) provenienti da osservazioni sperimentali, definiamo le differenze tra i valori rilevati \(y_i\) e i valori che, in corrispondenza delle \(x_i\), fornirebbe la funzione modello \eqref{eq:50}. Per individuare queste differenze che chiameremo comunque residui abbiamo diverse possibilità (si vedano, per es., le scelte fatte nel modello geometrico e in quello algebrico nella sezione riguardante la regressione circolare) tra le quali, quella seguita finora, pone i residui uguali alla distanza verticale tra il dato \(y_i\) e la corrispondente previsione del modello. Di conseguenza le grandezze incognite dipendenti dai parametri \(a,\,b,\,\dots\) \begin{equation} e_i=y_i-f(x_i,a,b,\dots) \qquad i=1,\,2,\,\dots n\tag{51}\label{eq:51} \end{equation} rappresentano, appunto, le differenze del valore sperimentale \(y_i\) con il valore "previsto" dalla funzione modello \(f(x_i,a,b,\dots)\).
  3. È ora possibile costruire la funzione obiettivo ottenuta sommando i quadrati di tutti i residui \begin{equation} s^2(a,b,\dots)=\sum_{i\,=\,1}^n\Bigl[y_i-f(x_i,a,b,\dots)\Bigr]^2.\tag{52}\label{eq:52} \end{equation} Questa funzione dipende dagli stessi parametri, ancora incogniti, che caratterizzano i residui e, per come è stata definita, non può che avere segno positivo o nullo.
  4. Coerentemente con il criterio fondamentale già enunciato, ricerchiamo il minimo di questa funzione per mezzo di metodi analitici quando, poste uguali a zero le derivate parziali prime, risulta risolvibile il sistema che ne deriva. In alternativa, e più spesso, dovremo utilizzare metodi numerici che invece fanno uso di opportuni algoritmi. Sia l'approccio analitico che quello numerico restituiscono comunque quei valori dei parametri \(\hat a,\,\hat b,\,\dots\) detti parametri di best-fit, che garantiscono (trascuriamo il fatto che la funzione obiettivo possa presentare più minimi) il raggiungimento del minimo della funzione obiettivo: in termini formali, i parametri di best-fit devono perciò soddisfare la condizione \(s^2(\hat{a},\hat b,\dots)<s^2(a,b,\dots),\, \forall\, (a\not=\hat a\wedge b\not=\hat b\wedge\dots)\).
  5. La sostituzione nella funzione di fit \eqref{eq:50} dei parametri \(\hat a,\,\hat b,\,\dots\) ci permette di ottenere quella particolare funzione detta, appunto, funzione di best-fit e quindi, sovrapponendola alle osservazioni, poterne verificare visualmente l'adeguatezza oltreché, naturalmente, interpretarne fisicamente i parametri.
  6. Infine, la conoscenza della funzione di best-fit, permette il calcolo effettivo dei residui \begin{equation} e_i=y_i-f(x_i,\hat a,\hat b,\dots) \qquad i=1,\,2,\,\dots n\tag{53}\label{eq:53} \end{equation} la cui distribuzione, come visto, può fornire non solo utili informazioni sul modello scelto ma pure, attraverso la stima della loro deviazione standard \(\sigma\),
  7. permettere il calcolo, analitico o numerico, delle deviazioni standard dei parametri \begin{equation} \sigma_{\hat a},\,\sigma_{\hat b},\dots\ .\tag{54}\label{eq:54} \end{equation}

Minimi quadrati e funzioni non linearizzabili

In questa sezione intendiamo presentare brevemente due esempi di applicabilità del metodo dei minimi quadrati nel caso di funzioni modello non lineari e, diversamente da quanto riportato nella sezione sui modelli alternativi, non linearizzabili. Per il calcolo dei risultati ci affidiamo unicamente a procedure numeriche per cui evitiamo di proporre formule risolutive peraltro inesistenti nel secondo esempio.

Le distribuzioni delle misure da cui partiamo sono riportate in figura 26. Questa, nella parte sinistra, suggerisce un andamento parabolico mentre le misure della parte destra mostrano un evidente andamento periodico oscillatorio.

Fig. 26. Modelli non linearizzabili tra le grandezze \(X\) e \(Y\).

Nel primo caso, scegliamo come funzione di fit del modello l'equazione quadratica \begin{equation} y=f(x,a,b,c)=ax^2+bx+c \end{equation} dipendente da tre parametri \(a,\,b,\,c\). Di conseguenza, come applicazione della \eqref{eq:52}, la funzione obiettivo diviene \begin{equation} s^2(a,b,c)=\sum_{i\,=\,1}^n\Bigl[y_i-(ax_i^2+bx_i+c)\Bigr]^2=\sum_{i\,=\,1}^n\Bigl(y_i-ax_i^2-bx_i-c\Bigr)^2\tag{55}\label{eq:55} \end{equation} anch'essa dipendente dagli stessi tre parametri.

Allo stesso modo e per i dati riportati nella parte destra di fig. 26, ipotizziamo una funzione di fit espressa da una combinazione lineare delle funzioni goniometriche \(a\sin(\omega\, x)\) e \(b\cos(\omega\, x)\) \begin{equation} y=f(x,a,b,c,\omega)=a \sin(\omega\, x)+b \cos(\omega\, x)+c. \end{equation} In tal caso sono quattro i parametri incogniti. La funzione obiettivo collegata alle misure e al modello scelto assume la forma \begin{equation} s^2(a,b,c,\omega)=\sum_{i\,=\,1}^n\Bigl[y_i-a \sin(\omega\, x_i)-b \cos(\omega\, x_i)-c\Bigr]^2. \end{equation} In entrambi i modelli, il calcolo numerico condotto con strumenti presenti anche in un normale foglio elettronico fornisce i valori dei parametri ottimali e le relative incertezze. Di conseguenza diviene possibile tracciare le funzioni di best-fit \[ y=\hat a x^2+\hat b x+\hat c,\hskip2cm y=\hat a\sin(\hat \omega\,x)+\hat b\cos(\hat \omega\,x)+\hat c \] sovrapponendole alla distribuzione delle misure (fig. 26).
Considerando la presenza di \(p=3\) parametri nel primo esempio e \(p=4\) nel secondo e in base alla \[ \sigma^2={1\over n-p}\sum_{i\,=\,1}^n e_i^2 ={1\over n-p}\sum_{i\,=\,1}^n\Bigl[y_i-f(x_i,\hat a,\hat b,\dots)\Bigr]^2, \] evidente generalizzazione della \eqref{eq:24}, associamo infine alla distribuzione dei residui pure la stima della loro deviazione standard (fig. 27).

Fig. 27. Distribuzione dei residui e loro deviazione standard.

Regressione circolare

In questa sezione intendiamo affrontare ancora un caso particolare di modello non lineare dove la funzione scelta sia rappresentativa di una circonferenza. In altre parole vorremmo determinare l'equazione di quella circonferenza che si adatti nel modo migliore ad un insieme di osservazioni sperimentali rappresentate, al solito, dalle coppie di misure \((x_i,y_i)\) con \(i=1,2,\dots,n\).

Come riportato nella pagina associata a questa un tale problema, pur meno frequente del modello lineare, si incontra in diverse occasioni, anche in un laboratorio scolastico. Per esempio può sorgere la necessità di misurare, sulla base di fotografie, il diametro di un massimo di diffrazione da foro circolare (diffrazione) oppure il raggio della traiettoria circolare di un fascio di elettroni in moto in un campo magnetico (traiettoria circolare).

Sulla base dello stesso insieme di dati e della medesima funzione modello, l'equazione cartesiana di una circonferenza, proponiamo due modi alternativi per definire i residui: indichiamo il primo approccio come modello geometrico e il secondo come modello algebrico. Per entrambi costruiamo la conseguente funzione obiettivo \(s^2\) e, al variare dei parametri dai quali questa dipende, ricerchiamo le condizioni di minimo e ne confrontiamo i risultati.

Fig. 28. Circonferenza ottimale.

Modello geometrico. Anche con riferimento alla figura 28 dove riportiamo assieme alla distribuzione delle misure pure una circonferenza (è quella ottimale), sia \(\gamma\) l'equazione che descrive analiticamente la circonferenza di centro \(C=(a,b)\) e raggio \(r\) e che rappresenta pure la funzione modello, \[ \gamma: (x-a)^2+(y-b)^2=r^2.\tag{56}\label{eq:56} \] Ciascuna coppia sperimentale \((x_i,y_i)\) possiede dal centro \(C\) una distanza \(d_i\) pari a \[ d_i=\sqrt{(x_i-a)^2+(y_i-b)^2} \] mentre la distanza \(e_i\) del medesimo punto dalla circonferenza è invece espressa dalla differenza \begin{equation} e_i=d_i-r=\sqrt{(x_i-a)^2+(y_i-b)^2}-r\tag{57}\label{eq:57} \end{equation} Quest'ultima grandezza possiede una interpretazione geometrica immediata in quanto appare come la lunghezza del segmento di direzione radiale che collega la coppia \((x_i,y_i\)) alla circonferenza (fig. 28) e questa caratteristica dà ragione del nome di questo approccio.
Evidentemente queste differenze assumono segno positivo per i punti esterni alla circonferenza mentre il segno è invece negativo per quelli interni e, in questo contesto, appaiono l'equivalente dei residui definiti come distanze verticali tra la misura \(i\)-esima e la previsione del modello nello schema di figura 25.
La funzione obiettivo associata a questa scelta è quindi data dalla somma dei quadrati di questi residui \begin{equation} s^2(a,b,r)=\sum_{i\,=\,1}^n e_i^2=\sum_{i\,=\,1}^n \Bigl[\sqrt{(x_i-a)^2+(y_i-b)^2}-r\Bigr]^2\tag{58}\label{eq:58} \end{equation} e mostra di dipendere dai tre parametri \(a,\,b,\,r\). Per dedurre i valori ottimali \(\hat a,\,\hat b,\,\hat r\) in corrispondenza dei quali \(s^2(a,b,r)\) presenta un minimo, sappiamo che ciò comporta risolvere analiticamente il sistema ottenuto annullando le tre derivate parziali della \eqref{eq:58}. La presenza della radice rende comunque questo calcolo impraticabile per cui dobbiamo affidarci ad algoritmi di carattere numerico: questi, nel caso le misure si distribuiscano lungo tutta la circonferenza forniscono risultati sicuramente affidabili e, nel nostro esempio, riportati nella fig. 28. Nella figura 29 proponiamo invece la distribuzione dei residui con la relativa deviazione standard: l'indice posto in ascissa assume il suo valore iniziale in corrispondenza del dato che in figura 28 forma, con l'asse passante per il centro e parallelo all'asse \(X\), l'angolo più prossimo allo zero.

Fig. 29. Residui della circonferenza ottimale di figura 28.

Modello algebrico. La presenza di radicali nella funzione obiettivo \eqref{eq:58} comporta indubbi problemi di carattere analitico per cui cerchiamo per tale funzione una definizione alternativa. A tale scopo notiamo che se una coppia di dati \((x_i,y_i)\) appartiene alla circonferenza \(\gamma\) allora essa soddisfa ovviamente l'equazione \[ (x_i-a)^2+(y_i-b)^2-r^2=0, \] mentre se la sua distanza dal centro è maggiore del raggio \(r\), allora dev'essere \[ (x_i-a)^2+(y_i-b)^2-r^2>0. \] Nel caso opposto, quanto la distanza di \((x_i,y_i)\) da \(C\) è minore di \(r\) risulta invece \[ (x_i-a)^2+(y_i-b)^2-r^2<0. \] Abbiamo quindi a disposizione una grandezza che si comporta come i residui definiti nel modello precedente \eqref{eq:57}. Seppure questa quantità sia difficilmente interpretabile geometricamente (da ciò l'aggettivo "algebrico" per questo approccio) possiamo sostituire la definizione dei residui con la seguente \begin{equation} e_i=(x_i-a)^2+(y_i-b)^2-r^2\tag{59}\label{eq:59} \end{equation} per cui la funzione obiettivo, in questa ipotesi, diviene \begin{equation} s^2(a,b,r)=\sum_{i\,=\,1}^n e_i^2=\sum_{i\,=\,1}^n \Bigl[{(x_i-a)^2+(y_i-b)^2}-r^2\Bigr]^2.\tag{60}\label{eq:60} \end{equation}

Quest'ultima funzione non presenta radicali e possiede il non trascurabile vantaggio che sviluppando i calcoli delle tre derivate e con il cambio di variabili,

\begin{equation} z_i=x_i^2+y_i^2,\qquad B=-2a,\qquad C=-2b,\qquad D=a^2+b^2-r^2,\tag{61}\label{eq:61} \end{equation}

il sistema che ne deriva è lineare e quindi risolubile analiticamente (applichiamo questa "variante" al modello algebrico nella "demonstration" alla pagina successiva). Il calcolo numerico condotto sulla funzione obiettivo \eqref{eq:60} oppure sul sistema lineare fornisce risultati che, a meno di un centesimo, coincidono con quanto riportato nelle precedenti due figure.
Per un confronto più preciso tra gli esiti di entrambi i modelli proponiamo infine le coppie di figure seguenti. Queste, a partire dallo stesso insieme di misure, confermano come le differenze per centro e raggio siano dell'ordine dei centesimi. Pertanto i due modelli sono entro questi limiti equivalenti, sebbene l'odierna diffusione del calcolo numerico sembra preferire l'approccio geometrico su quello algebrico.

Fig. 30. Regressione circolare: confronto tra i modelli geometrico e algebrico.
Fig. 31. Regressione circolare: modelli geometrico e algebrico.

Bibliografia

[1] N. Chernov, Circular and Linear Regression, CRC Press, (2011)

[2] S. Chatterjee, A. Hadi, Regression Analysis by Example, Wiley-Interscience, (2006)

[3] D. Montgomery, E. Peck, G. Vining, Linear Regression Analysis, Wiley, (2012)

[4] J. Taylor, Introduzione all'analisi degli errori, Zanichelli, (1986)

[5] P. Bevington, Data Reduction and Error Analysis for the Physical Sciences, McGraw-Hill, (1969)

Segue la pagina sull'utilizzo didattico della regressione