La formula per il fit di un set di dati con una retta usando SymPy

In [3]:
from sympy import *

La somma degli scarti di dati $(x_i,y_i)$ (pesati con l'errore $\sigma_i$ sugli $y_i$) rispetto una linea $y = b \cdot x + a$ è: $$ \chi^2 = \sum_i \frac{(y_i - b \cdot x_i - a)^2}{\sigma_i^2} $$

Per trovare il minimo di questa espressione al variare di $b$ e $a$, cerchiamo gli zeri contemporanei delle derivate di $\chi^2$ rispetto ad $a$ e rispetto a $b$, generalizzando in modo ovvio la formula per gli estremi in una variabile.</br> È sufficiente limitarsi al caso di tre sole misure e poi generalizzare.

Domanda: Perchè tre misure e non una oppure due?

In [16]:
x1,x2,x3,y1,y2,y3,a,b,sq1,sq2,sq3 = symbols('x1 x2 x3 y1 y2 y3 a b sq1 sq2 sq3')
In [17]:
chi2 = (y1 - b*x1 - a)**2/sq1 + (y2 - b*x2 - a)**2/sq2 + (y3 - b*x3 - a)**2/sq3

Derivata rispetto ad a:

In [18]:
chi2_prime_a = diff(chi2,a)
chi2_prime_a
Out[18]:
$\displaystyle \frac{2 a + 2 b x_{3} - 2 y_{3}}{sq_{3}} + \frac{2 a + 2 b x_{2} - 2 y_{2}}{sq_{2}} + \frac{2 a + 2 b x_{1} - 2 y_{1}}{sq_{1}}$

Espandiamo il risultato per poi estrarre i coefficienti di a e b:

In [26]:
chi2_prime_a_expand = expand(chi2_prime_a)
chi2_prime_a_expand
Out[26]:
$\displaystyle \frac{2 a}{sq_{3}} + \frac{2 a}{sq_{2}} + \frac{2 a}{sq_{1}} + \frac{2 b x_{3}}{sq_{3}} + \frac{2 b x_{2}}{sq_{2}} + \frac{2 b x_{1}}{sq_{1}} - \frac{2 y_{3}}{sq_{3}} - \frac{2 y_{2}}{sq_{2}} - \frac{2 y_{1}}{sq_{1}}$
In [29]:
collect(chi2_prime_a_expand,(a,b))
Out[29]:
$\displaystyle a \left(\frac{2}{sq_{3}} + \frac{2}{sq_{2}} + \frac{2}{sq_{1}}\right) + b \left(\frac{2 x_{3}}{sq_{3}} + \frac{2 x_{2}}{sq_{2}} + \frac{2 x_{1}}{sq_{1}}\right) - \frac{2 y_{3}}{sq_{3}} - \frac{2 y_{2}}{sq_{2}} - \frac{2 y_{1}}{sq_{1}}$

Derivata rispetto a b:

In [19]:
chi2_prime_b = diff(chi2,b)
chi2_prime_b
Out[19]:
$\displaystyle - \frac{2 x_{3} \left(- a - b x_{3} + y_{3}\right)}{sq_{3}} - \frac{2 x_{2} \left(- a - b x_{2} + y_{2}\right)}{sq_{2}} - \frac{2 x_{1} \left(- a - b x_{1} + y_{1}\right)}{sq_{1}}$

Espandiamo e estraiamo i coefficienti di a e b:

In [30]:
chi2_prime_b_expand = expand(chi2_prime_b)
chi2_prime_b_expand
Out[30]:
$\displaystyle \frac{2 a x_{3}}{sq_{3}} + \frac{2 a x_{2}}{sq_{2}} + \frac{2 a x_{1}}{sq_{1}} + \frac{2 b x_{3}^{2}}{sq_{3}} + \frac{2 b x_{2}^{2}}{sq_{2}} + \frac{2 b x_{1}^{2}}{sq_{1}} - \frac{2 x_{3} y_{3}}{sq_{3}} - \frac{2 x_{2} y_{2}}{sq_{2}} - \frac{2 x_{1} y_{1}}{sq_{1}}$
In [31]:
collect(chi2_prime_b_expand,(a,b))
Out[31]:
$\displaystyle a \left(\frac{2 x_{3}}{sq_{3}} + \frac{2 x_{2}}{sq_{2}} + \frac{2 x_{1}}{sq_{1}}\right) + b \left(\frac{2 x_{3}^{2}}{sq_{3}} + \frac{2 x_{2}^{2}}{sq_{2}} + \frac{2 x_{1}^{2}}{sq_{1}}\right) - \frac{2 x_{3} y_{3}}{sq_{3}} - \frac{2 x_{2} y_{2}}{sq_{2}} - \frac{2 x_{1} y_{1}}{sq_{1}}$

A questo punto è facile scrivere l'equazione matriciale che determina $a$ e $b$:

$$ \begin{bmatrix} \sum \frac{1}{\sigma_i^2} & \sum \frac{x_i}{\sigma_i^2} \\ \sum \frac{x_i}{\sigma_i^2} & \sum \frac{x^2_i}{\sigma_i^2} \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum \frac{y_i}{\sigma_i^2} \\ \sum \frac{x_i y_i}{\sigma_i^2} \end{bmatrix} $$$$ \begin{bmatrix} Sinv & Sx \\ Sx & Sxsq \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} Sy \\ Sxy \end{bmatrix} $$
In [6]:
a, b, Sinv, Sx, Sxsq, Sy, Sxy = symbols('a b Sinv Sx Sxsq Sy Sxy')

eq1 = Sinv*a + Sx*b - Sy
eq2 = Sx*a  + Sxsq*b -Sxy

solve((eq1,eq2),a,b)
Out[6]:
{a: (-Sx*Sxy + Sxsq*Sy)/(Sinv*Sxsq - Sx**2),
 b: (Sinv*Sxy - Sx*Sy)/(Sinv*Sxsq - Sx**2)}

Nel caso particolare in cui tutti i $\sigma_i$ sono uguali, definendo $\bar{u} = \sum u_i /N$ ($u = 1,x,y,xy$) l'equazione che determina $a, b$ diventa:

$$ M \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} 1 & \overline{x_i} \\ \overline{x} & \overline{x^2} \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \overline{y} \\ \overline{xy} \end{bmatrix}, \qquad M^{-1} = \frac{1}{\Delta} \begin{bmatrix} \overline{x^2} & -\overline{x_i} \\ -\overline{x} & 1 \end{bmatrix},\,\, \Delta=\overline{x^2}-\overline{x}^2 $$


con soluzione $a = (\overline{x^2}\cdot \overline{y}-\overline{x}\cdot\overline{xy})/\Delta$, $b = (\overline{xy}-\overline{x}\cdot\overline{y})/\Delta$


È facile mostrare che il punto $(\overline{x},\overline{y})$ appartiene alla retta di fit cioè che $ \overline{y} = b \cdot \overline{x} +a $.

In [ ]: