CARAMELLE di Forcone (ψ)

Supponiamo di avere una serie di osservazioni $(x_i, y_i \pm \sigma_i)$ e di voler trovare la “migliore” retta $y = a x + b$ che spieghi queste osservazioni. Questo problema è noto come adattamento (in inglese fit). Chiaramente non tutti i dati sono spiegabili da una retta, ma molto spesso ci si può ridurre a questo caso anche se i dati non sono lineari, perciò questo caso è particolarmente importante. Ad esempio se abbiamo una relazione esponenziale $y = a e^{bx}$ prendendo il logaritmo da entrambi i lati abbiamo $\log{y} = \log{a} + b x$ e quindi le osservazioni $(x_i, \log{y_i})$ seguono una linea retta.

Per ottenere una risposta univoca bisogna decidere cosa intendiamo per retta “migliore”. Per motivi statistici che non andiamo ad indagare in questo caso per “miglior retta” intendiamo quella che minimizza la quantità $\chi^2$ definita come

$$\chi^2 \equiv \sum_{i=1}^N \frac{\pqty{f(x_i)-y_i}^2}{\sigma_i^2}$$

dove $f(x) = ax+b$, in questo caso una retta, è la funzione che riteniamo possa ben rappresentare i dati. Innanzitutto notiamo che $\chi^2$ è adimensionale; inoltre se $f$ approssima i dati molto bene il numeratore sarà piccolo e quindi $\chi^2$ sarà piccolo, perciò vogliamo minimizzarla. La varianza nel denominatore ci dà una misura di controllo per sapere se $f(x_i)-y_i$ è da considerarsi grande o piccolo: in particolare i punti poco precisi (che hanno un errore $\sigma_i$ grande) conteranno poco nella somma, mentre quelli più precisi conteranno di più.

Nel caso della retta siamo in grado di effettuare la minimizzazione in maniera analitica. Consideriamo $\chi^2$ come una funzione di $a$ e $b$ e la minimizziamo calcolandone le derivate. Abbiamo

$$\begin{align} \pdv{\chi^2}{a} &= \sum_{i=1}^N \frac{2 x_i \pqty{ax_i+b-y_i}}{\sigma_i^2}=0\\ \pdv{\chi^2}{b} &= \sum_{i=1}^N \frac{2\pqty{ax_i+b-y_i}}{\sigma_i^2}=0 \end{align}$$

Esplicitando $a$ e $b$ le due equazioni sono

$$\begin{align} a\sum_{i=1}^N \frac{x_i^2}{\sigma_i^2} + b\sum_{i=1}^N \frac{x_i}{\sigma_i^2} -\sum_{i=1}^N \frac{x_i y_i}{\sigma_i^2}=0\\ a\sum_{i=1}^N \frac{x_i}{\sigma_i^2} + b\sum_{i=1}^N \frac{1}{\sigma_i^2} -\sum_{i=1}^N \frac{y_i}{\sigma_i^2}=0 \end{align}$$

I risultati sono ottenuti nella maniera più concisa definendo l’operazione

$$K [U] \equiv \sum_{i=1}^N \frac{U_i}{\sigma_i^2}$$

per una quantità generica $U$. Con questa notazione, $\chi^2 = K[(y-ax-b)^2]$. Perciò troviamo

$$\begin{align} a &= \frac{K[xy] K[1] -K[x] K[y]}{K[x^2] K[1]-K[x]^2} \\ b &= \frac{K[x^2] K[y] -K[xy] K[x]}{K[x^2] K[1]-K[x]^2} \end{align}$$

Una volta ottenute queste formule possiamo anche calcolare l’errore su $a$ e $b$. Infatti vicino al minimo di una funzione di due variabili $z = (a,b)$ avremo

$$\chi^2(z) \approx \chi^2(z^*) + \frac12 \sum_{i,j=1}^2 (z_i -z_i^*) \pdv{^2 \chi^2}{z_i \partial z_j} (a^*, b^*) (z_j -z_j^*)$$

dove $H$ è l’Hessiana di $\chi^2$ (cioè la matrice delle derivate parziali) e $(a^*, b^*)$ sono i valori ottimali che abbiamo trovato sopra. La matrice di covarianza di $a$ e $b$ sarà perciò $\Sigma$ tale che

$$(\Sigma^{-1})_{ij} =\frac12 \pdv{^2 \chi^2}{z_i \partial z_j} (a^*, b^*)$$

Il calcolo è lungo e noioso ma elementare e troviamo che

$$\Sigma = \begin{pmatrix} \sigma^2(a) & \mathrm{Cov}(a,b)\\ \mathrm{Cov}(a,b) & \sigma^2(b) \end{pmatrix}$$

con le seguenti varianze/covarianze:

$$\begin{align} \sigma^2(a) &= \frac{K[1]}{K[x^2] K[1]-K[x]^2}\\ \sigma^2(b) &= \frac{K[x^2]}{K[x^2] K[1]-K[x]^2}\\ \mathrm{Cov}(a,b) &= -\frac{K[x]}{K[x^2] K[1]-K[x]^2} \end{align}$$

È curioso notare che nel caso in cui la retta $y=ax+b$ degenera ad una costante, ovvero se poniamo $a=0$ e cerchiamo la “migliore costante”, un semplice calcolo mostra che il valore di $b$ meglio adattato è

$$b = \frac{\sum_{i=1}^N y_i / \sigma_i^2}{\sum_{i=1}^N 1 / \sigma_i^2}$$

che esattamente la media ponderata degli $y$ che minimizza l’errore sulla media, come abbiamo visto in un precedente articolo.