Queste note sono basate sugli appunti fatti con
Gianmarco Manzini negli anni 1995-2005
Metodo QR per sistemi lineari
Appunti di calcolo numerico
Autore/Autrice
Affiliazione
Enrico Bertolazzi
University of Trento, Department of Industrial Engineering
Problema dei minimi quadrati
Consideriamo il seguente sistema lineare
\left\{
\begin{aligned}
x &= 0 \\
x &= 1
\end{aligned}
\right.
È evidente che questo sistema non ha soluzione. Supponiamo ora che x rappresenti una quantità ottenuta tramite misurazioni soggette a errore, e che x = 0 e x = 1 siano i valori di due misure. Possiamo quindi riformulare il problema come
\left\{
\begin{aligned}
x &= 0+\epsilon_1 \\
x &= 1+\epsilon_2
\end{aligned}
\right.
dove \epsilon_1 ed \epsilon_2 rappresentano gli errori di misura. Poiché tali errori non sono noti, possiamo distribuire l’errore tra le due misure, ottenendo
\left\{
\begin{aligned}
x &= 0+\epsilon_1 \\
x &= 1+\epsilon_2 \\
x &= 1+\epsilon_3
\end{aligned}
\qquad
\color{red}
\right\}
\quad
\textrm{\color{blue}come distribuisco gli errori?}
alcune possibilità
minimizzo |\epsilon_1|+|\epsilon_2|+|\epsilon_3|, ma la funzione valore assoluto non è differenziabile sullo 0.
minimizzo \epsilon_1^2+\epsilon_2^2+\epsilon_3^2, questa funzione è differenziabile.
posso migliorare se minimizzo w_1\epsilon_1^2+w_2\epsilon_2^2+w_3\epsilon_3^2, dove w_k sono dei pesi positivi che stanno a indicare la affidabilità della misuira. Se w_k è grande la misura è affidabile se w_k è piccolo lo è meno.
g'(x) = 2x+4(x-1) = 6x-4 = 0,\quad\implies\quad
x = \dfrac{2}{3}
Generalizzando a n misure:
\left\{
\begin{aligned}
x &= y_1+\epsilon_1 \\
x &= y_2+\epsilon_2 \\
x &= y_3+\epsilon_3 \\
&\vdots \\
x &= y_n+\epsilon_n
\end{aligned}
\right.
che è la media delle misure. Soluzione molto intuitiva.
Questo approccio si può generalizzare e prende il nome di metodo dei minimi quadrati.
Regressione lineare
Consideriamo una serie di punti (x_k,y_k) e supponiamo di sapere che questi punti si discribuiscono su una retta
R(x) = A + Bx
ovvero y_k = R(x_k) con A e B da determinare. Questa si chiama retta di regressione lineare. Nel caso di dati perfetti basterebbero due punti ma i dati sono affetti da errore e quindi
\left\{
\begin{aligned}
A + B x_0 &= y_0 + \epsilon_0 \\
A + B x_1 &= y_1 + \epsilon_1 \\
A + B x_2 &= y_2 + \epsilon_2 \\
&\vdots \\
A + B x_n &= y_n + \epsilon_n
\end{aligned}
\right.
per cui la somma dei quadrati degli errori diventa
\sum_{k=0}^n \epsilon_k^2 =
\sum_{k=0}^n (A + B x_k - y_k)^2 =
f(A,B)
i punti di minimo soddisfano \nabla f(A,B)=0 ovvero
\begin{aligned}
\dfrac{\partial f(A,B)}{\partial A} &=
2\sum_{k=0}^n (A + B x_k - y_k) \\
& =
2(n+1)A + 2\left(\sum_{k=0}^n x_k\right) B -
2\sum_{k=0}^n y_k = 0
\\
\dfrac{\partial f(A,B)}{\partial B} &=
2\sum_{k=0}^n (A + B x_k - y_k)x_k \\
& =
2\left(\sum_{k=0}^n x_k\right) A+
2\left(\sum_{k=0}^n x_k^2\right) B -
2\sum_{k=0}^n x_ky_k = 0
\end{aligned}
la cui soluzione produce la retta di regressione lineare cercata.
Usando la libreria di R per fare il fitting otteniamo
# Generazione di dati sinteticiset.seed(123) # Per riproducibilitàx <-seq(1, 10, by =0.5) # Variabile indipendentey <-2.5* x +rnorm(length(x), mean =0, sd =3) # Variabile dipendente con rumore# Fitting linearefit <-lm(y ~ x)# Sommario del modellosummary(fit)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-6.2312 -1.9646 -0.0946 1.4675 4.9818
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.7839 1.5685 0.500 0.624
x 2.4524 0.2553 9.606 2.78e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.047 on 17 degrees of freedom
Multiple R-squared: 0.8444, Adjusted R-squared: 0.8353
F-statistic: 92.28 on 1 and 17 DF, p-value: 2.781e-08
e poi plottiamo la soluzione
# Plot dei dati e della retta di regressioneplot(x, y, main ="Regressione lineare",xlab ="X", ylab ="Y", pch =16, col ="blue")abline(fit, col ="red", lwd =2)# Aggiunge la retta di regressione# Aggiunta della legendalegend("topleft", legend =c("Dati osservati", "Fit lineare"),col =c("blue", "red"), pch =c(16, NA), lwd =c(NA, 2), bty ="n")
Consideriamo una serie di punti (x_k,y_k) e supponiamo di sapere che questi punti si discribuiscono su una parabola
R(x) = A + Bx + C x^2
ovvero y_k = R(x_k) con A, B e C da determinare. Nel caso di dati perfetti basterebbero tre punti ma i dati sono affetti da errore e quindi
\left\{
\begin{aligned}
A + B x_0 + C x_0^2 &= y_0 + \epsilon_0 \\
A + B x_1 + C x_1^2 &= y_1 + \epsilon_1 \\
A + B x_2 + C x_2^2 &= y_2 + \epsilon_2 \\
&\vdots \\
A + B x_n + C x_n^2 &= y_n + \epsilon_n
\end{aligned}
\right.
per cui la somma dei quadrati degli errori diventa
\dfrac{1}{2}\sum_{k=0}^n \epsilon_k^2 =
\dfrac{1}{2}\sum_{k=0}^n (A + B x_k + C x_k^2 - y_k)^2 =
f(A,B,C)
i punti di minimo soddisfano \nabla f(A,B)=0 ovvero
\begin{aligned}
\dfrac{\partial f(A,B,C)}{\partial A} &=
\sum (A + B x_k + C x_k^2 - y_k) \\
& =
(n+1)A +
\left(\sum x_k\right) B +
\left(\sum x_k^2\right) C -
\sum y_k = 0
\\
\dfrac{\partial f(A,B,C)}{\partial B} &=
\sum (A + B x_k + C x_k^ 2- y_k)x_k \\
& =
\left(\sum x_k\right) A+
\left(\sum x_k^2\right) B +
\left(\sum x_k^3\right) C -
\sum x_ky_k = 0
\\
\dfrac{\partial f(A,B,C)}{\partial C} &=
\sum (A + B x_k + C x_k^ 2- y_k)x_k \\
& =
\left(\sum x_k^2\right) A+
\left(\sum x_k^3\right) B +
\left(\sum x_k^4\right) C -
\sum x_k^2y_k = 0
\end{aligned}
