%% 46   Minsta kvadratmetoden och kurvanpassning
% Innehåll
% 46.1   Minsta kvadratmetoden
% 46.2   Anpassning av modell till data
%% 46.1   Minsta kvadratmetoden
%
% Ett poplulärt sätt att lösa ekvationen A\mathbf x=\mathbf y\,, där A\, och \mathbf y är kända, i Matlab är med kommandot "vänsterdivision" med A\,:
% x = A \ y
% Tanken är att det ska motsvara multiplikation med A^{-1}\, från vänster,  men A\, behöver inte vara inverterbar och inte ens kvadratisk för att A\mathbf x=\mathbf y ska vara lösbar. Ett fenomen är att  x = A \ y  ger ett resultat även när lösning saknas och bara étt svar även om det finns flera.
% Ta, för att begripa vad som händer, en enkvation A\mathbf x=\mathbf y\, som saknar lösning. Det betyder att \mathbf y inte ligger i kolonnrummet till A\, dvs \mathbf y\not\in\mbox{Col}\,A\,.
%
% Ortogonala projektionen \hat \mathbf y\, på \mbox{Col}\, A\,  ligger däremot i \mbox{Col}\,A\, och är den vektor i kolonnrummet som ligger närmast \mathbf y. Ekvationen A\mathbf x=\hat\mathbf y\, är den ekvation som Matlab löser med  x = A \ y . Finns det flera lösningar \mathbf x\, levereras en av de lösnigar vars koordinater innehåller flest nollor.
% Definiton   En lösning \hat \mathbf x\, till A\mathbf x=\hat \mathbf y där \hat\mathbf y är ortogonala projektionen av \mathbf y\, på \mbox{Col}\,A kallas en minsta kvadratlösning till ekvationen A\mathbf x=\mathbf y\,. Avståndet |\!|\mathbf y-\hat\mathbf y|\!|\, kallas minsta kvadratfelet.
%
% Kommandot
% x = A \ y
% levererar alltså en minsta kvadratlösning till  ekvationen A\mathbf x=\mathbf y\,. När ekvationen är lösbar är \hat\mathbf x\, en lösning till A\mathbf x=\mathbf y\, men inte annars.
% Det här kan du utnyttja för att bestämma vilken vektor i ett delrum W\, som ligger närmast en given vektor \mathbf v\, och bestämma avståndet mellan \mathbf v och W\,.
%
% Uppgift 46.1 Bestäm den vektor som ligger närmast
% \mathbf v=
% \left(\begin{array}{r}
% 13 \\ -31 \\ 15 \\ -34
% \end{array}\right)\,.
% i det  linjära höljet till vektorerna
% \mathbf w_1 =
% \left(
%      \begin{array}{r}
%      -16\\    7\\   -6\\   -8
%      \end{array}
%      \right)\,,\quad
% \mathbf w_2 =
% \left(
%      \begin{array}{r}
%      15\\  -11\\   32\\  -13
%      \end{array}
%      \right)\,,\quad
% \mathbf w_3 =
% \left(
%      \begin{array}{r}
%      3\\   15\\    0\\   -3
%      \end{array}
%      \right)\,.
% Bestäm också avståndet mellan \mathbf v\, och det linjära höljet.
% Lösning



%
% Men - tänker du - hur ska jag kunna klara detta med handräkning - jag vet ju inte hur Matlab gör för att hitta en minsta kvadratlösning till A\mathbf x=\mathbf y\,. Måste jag då gå hela vägen med att först bestämma en ortogonalbas för \mbox{Col}\,A\, sen ortogonal projektionen \hat \mathbf y\, av \mathbf y\, på \mbox{Col}\,A\,osv?
% Nej då - och detta är väldigt elegant - det finns en genväg! Häng med nu...
% För att hitta en minsta kvadrat lösning \hat\mathbf x till A\mathbf x=\mathbf y ska du bestämma \hat\mathbf x så att A\hat\mathbf x=\hat\mathbf y\in \mbox{Col}\,A\, men det inträffar precis när \mathbf y-\hat \mathbf y=\mathbf y-A\mathbf \hat\mathbf x\, är vinkelrät mot varje kolonn i A\, och det kan du skriva
% \mathbf 0 = A^\top(\mathbf y-A\hat\mathbf x) = A^\top\mathbf y-A^\top A\hat\mathbf x\,.
%
% Det betyder att minsta kvadratlsöningarna \hat \mathbf x\,till A\mathbf x=\mathbf y är precis lösningarna till
% A^\top A\mathbf x = A^\top \mathbf y
% och den ekvationen är alltid lösbar!
% Definition   Ekvationen A^\top A\mathbf x = A^\top \mathbf y\,, där A\,och \mathbf y är kända, kallas normalekvationen till A\mathbf x =\mathbf y.
% Du har alltså
% Sats   Minsta kvadratlösningarna till A\mathbf x=\mathbf y\,, där A\,och \mathbf y är kända, är precis lösningarna till normalekvationen A^\top A\mathbf x = A^\top \mathbf y\,.
% Du har nu en ny metod att beräkna ortogonala projektionen av en vektor \mathbf v\in \mathbb R^n på ett delrum W\, till \mathbb R^n\, som inte går vägen via en ortogonalbas för W\,.
% Metod (för att bestämma ortogonala projektionen av \mathbf v på W\,)
% Bestäm vektorer \mathbf w_1,\,\mathbf w_2,\,\ldots,\,\mathbf w_q som spänner ut W\,.
% Sätt A = \left(\begin{array}{rrr}
% \mathbf w_1 &\mathbf w_2 &\ldots &\mathbf w_q
%  \end{array}\right)\,.
% Bestäm én lösning \hat\mathbf x\, till normalekvationen A^\top A\mathbf x= A^\top\mathbf v\,.
% Då gäller att \mbox{proj}_W(\mathbf v)=\hat\mathbf v= A\hat\mathbf x\,.
% Uppgift 46.1   Låt W vara nollrummet till matrisen B\, där
% B = \left(\begin{array}{rrrr} -5 & 5 & 34 & 2\\ -17 & 16 & 112 & 10\\ 21 & -18 & -132 & -18 \end{array}\right)\,.
% Använd normalekvationen för att bestäm den ortogonala projektionen av vektorn \mathbf v\,  på  W\, när
% \mathbf v=   \left(\begin{array}{r} -49\\ -68\\ 33\\ 31 \end{array}\right)\,.



%% 46.2   Anpassning av modell till data
%
% Säg att du gör observationer i ett experiment där en variabel y, beror av värdet på en annan som du kallar x\,. Du oberverar y\, när x\, har värdena x_1,\,x_2,\,\ldots,\, x_n\, och noterar då vädena y_1,\,y_2,\, \ldots ,\, y_n\, på y\,.  Du för väntar dig ett samband mellan x\, och y\, som är linjärt:  att y= ax+b\, och vill från mätdata bestämma vad a\, och b\,är. Du plottar punkterna (x_1,\,y_1),\, (x_2,\,y_2),\,\ldots,\,(x_n,\,y_n)\, och ser figuren
%
% Du ser att punkterna inte ligger längs en rät linje men misstänker att det beror på att du inte lyckats mäta exakt nog och möjligen på andra störningar i experimentet.
% För att bestämma a\, och b\,ska du lösa ekvationen
% \left(\begin{array}{rr} \mathbf x &  \mathbf 1\end{array}\right)\left(\begin{array}{r} a\\ b\end{array}\right) = \mathbf y\qquad (1)
% där
% \mathbf x = \left(\begin{array}{r} x_1 \\ x_2\\ \vdots \\ x_n\end{array}\right)\,,\quad \mathbf 1 = \left(\begin{array}{r} 1\\ 1\\ \vdots \\ 1\end{array}\right)\,,\quad
% \mathbf y = \left(\begin{array}{r} y_1 \\ y_2\\ \vdots \\ y_n\end{array}\right)\,
% men inser att den ekvationen kommer att sakna lösning eftersom dina mätdata inte ligger längs en linje. Du har ett ekvationssystem som är överbestämt . Du skulle kunna kasta bort alla punkter utom två så att systemtet blir lösbart - men i så fall vilka? Alla punkter kan vara värdefulla. Du får idén att istället försöka bestämma a\, och b\, så att vektorn till vänster i (1) kommer så nära \mathbf y\, som möjligt. Det betyder att du ska bestämma en minsta kvadratlösning till ekvationen. Med Matlab kan du göra det med kommandot
% v = ( x ones( n, 1) ) \ y
% där koordinaterna i \mathbf v\, är de värden på a\, och b\, som du söker.
% Testa detta nu med
% x = [ 49.8801, 18.9669, 36.0980, 34.1996, 25.4898, 15.6875, 11.0054, 26.8445, 17.3640, 39.0310, 24.8145, 43.6624, 39.3692, 32.8410, 17.0742, 48.2954, 20.6129, 46.9832, 18.9508, 24.9426 ]' ;
% y = [ 135.5154, 48.3017, 95.1003, 88.0494, 68.7012, 35.5368, 24.6224, 69.3722, 43.3655, 102.3095, 68.5492, 123.9928, 110.5660, 91.6542, 40.0557, 136.0589, 52.5936, 131.2502, 44.6031, 62.3139]' ;



%
% En mer allmän situation kan vara att du på teoretiska grunder vet, eller misstänker, att variabeln y\, är en linjärkombination av kända  funktioner som beror på variabeln x\,:
% y = c_1\,f_1(x)+c_2\,f_2(x)+\cdots + c_n\,f_n(x)
% Ekvationen kallas för en matematisk modell av experimentet och från matvärden (samplade data) hoppas du kunna bestämma vikterna (parametrarna)  c_1,\, c_2,\, \ldots,\, c_n\,.
% Du utför experiementet och får samplade väden \mathbf x_s på variabeln x\, i en kolonnvektor och motsvarande värden  \mathbf y_s\, för y\,. För att bestämma parameterarna ska du lösa ekvationsystemet
% \left(\begin{array}{rrrr} f_1(\mathbf x_s) & f_2(\mathbf x_s) & \ldots &f_n(\mathbf x_s)\end{array}\right)\left(\begin{array}{ r} c_1\\ c_2\\ \vdots \\ c_n\end{array}\right)=\mathbf y_s\,.
% Matrisen
% X=\left(\begin{array}{rrrr} f_1(\mathbf x_s) & f_2(\mathbf x_s) & \ldots &f_n(\mathbf x_s)\end{array}\right)\,
% kallas i detta sammanhang för designmatrisen och kolonnerna
% \mathbf \beta =\left(\begin{array}{ r} c_1\\ c_2\\ \vdots \\ c_n\end{array}\right),\quad \mathbf y_s\,
% för parameter- respektive obeservationsvektorn. Ekvationssystemet blir då
% X\mathbf \beta = \mathbf y_s
% och är normalt (helst) överbestämt.
% Uppgift 46.2    I ett försök försöker du bestämma paramterarna \beta_1,\, \beta_2,\, \beta_3 i den matematiska modellen
% y = \beta_1\,x^2+ \beta_2\,x + \beta_3\,.
% Du får mätdatat
% xs = [ 5.0958 , 4.4808 , 13.3567 , 16.8878 , 6.8892 , 15.6104 , 13.5066 , 0.1343 , 12.0434 , 7.7354 , 18.3198 , 0.0230 , 9.2490 , 8.4870 , 9.2183 , 15.4032 , 6.4494 , 15.6948 , 9.4271 , 0.7153 ]' ;
% ys = [ 51.6779, 50.6563, 59.2604, 39.8897, 59.5261, 49.6913, 57.6187, 17.2878, 61.8836, 64.1055, 29.4665, 16.2919, 63.4726, 63.0054, 63.0589, 50.7076, 59.4302, 48.9711, 64.5388, 22.4745 ]' ;
% Bestäm paramtrarna i modellen.



%
%
%