Föreläsning 3
Minsta kvadratmetoden, MKM
• Hur kan man anpassa en funktion till en serie mätdata?
• Normalekvationer, NE
• Residual och felkvadratsumma
• Vad menas med i minsta kvadrat mening?
• Teoretisk motivering
• Centrering
• Ickelinjär modell
Exempel
En rät linje, y = kx + m ska anpassas till följande mätvärden x 1 2 3
y 4 2 1
Hur kan man anpassa en funktion till en serie mätdata?
Ställ upp problemet på matrisform Ac = y
1 1 2 1 3 1
km
=
4 2 1
Detta är ett överbestämt system, dvs det är fler ekvationer än obekanta.
För att lösa ett sådant kan man förlänga med AT
1 2 3 1 1 1
1 1 2 1
km
=
1 2 3 1 1 1
4 2
⇒
Normalekvationer, NE
14 6 6 3
km
=
117
Det system som erhållits kallas normalekvationer och skrivs allmänt
ATAc = ATy c kan nu lösas ut
c =
−1.5000 5.3333
Plot av lösning
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
1 1.5 2 2.5 3 3.5 4
y
En rät linje anpassas till mätdata
Felet för första punkten
Felet för andra punkten
Felet för tredje punkten
* Mätdata
Residual och felkvadratsumma
Residualen är vektorn r = Ac − y
r = Ac − y ⇒
1 1 2 1 3 1
−1.5000 5.3333
−
4 2 1
=
−0.1667 0.3333
−0.1667
Felkvadratsumman
i=1 r2i kan beräknas med t.ex. rTr.
fks = (−0.1667)2 + 0.33332 + (−0.1667)2 = 0.1667
Lösningen som erhållits är den som minimerar residualvektorns euklidiska norm. Det innebär att felkvadratsumman inte kan bli mindre med ett annat val av k och m. Man säger att funktionen är anpassad till mätdata i minsta kvadratenmening.
Teoretisk motivering (ett alternativ till kursbokens)
En rät linje y = kx + m ska anpassas till mätdata (xi, yi), i = 1, 2, ..., n.
k och m ska bestämmas så att felkvadratsumman
fks(k, m) =
n i=1
(kxi + m − yi)2
minimieras.
Minimum för en funktion hittar man om man deriverar funktionen och sätter derivatorna lika med noll ...
Teoretisk motivering (ett alternativ till kursbokens)
∂f ks(k,m)
∂k = 2n
i=1 xi(kxi + m − yi)
∂f ks(k,m)
∂m = 2 n
i=1(kxi + m − yi) ⇒
∂f ks(k,m)
∂k = 0
∂f ks(k,m)
∂m = 0 ⇒
mn
i=1 xi + k n
i=1 x2i = n
i=1 xiyi mn + k n
i=1 xi = n
i=1 yi Detta är normalekvationerna!
Ytterligare exempel
Funktionen F (T ) = L0 + L1T där L1 = λL0 ska anpassas med MKM till följande mätvärden
Temp (C) 20.0 25.5 30.2 36.8 41.0 Längd (cm) 8.78 8.93 9.06 9.25 9.40
Bestäm λ.
Ytterligare exempel, matrisformulering och NE.
Ställ upp problemet på matrisform Ac = y
1 20.0 1 25.5 1 30.2 1 36.8 1 41.0
L0 L1
=
8.78 8.93 9.06 9.25 9.40
Formulera normalekvationerna ATAc = ATy
5.00 153.5 153.5 4997.53
L0 L1
=
45.42 1402.727
Ytterligare exempel, lösning
Lös ut c och beräkna λ !
c =
8.1866 0.0292
Vilket ger L0 = c1 = 8.1866 och λ = L1/L0 = c2/c1 = 0.0036
Ytterligare exempel, avrundar i NE.
Avrunda till tre korrekta siffror i normalekvationerna ATAc = ATy
5 154 154 5000
L0 L1
=
45.4 1400
Lös ut c och beräkna λ !
c =
8.8785 0.0065
Vilket ger L0 = c1 = 8.8785 och λ = L1/L0 = c2/c1 = 0.0007 Värdet på λ skiljer sig markant från första uträkningen!
Varför? Normalekvationerna är illakonditionerade.
Konditionstalet för ATA är 20000. Åtgärd: centrering!!
Ytterligare exempel, centrering
Mätvärdena anpassas istället till F (T ) = c1 + c2(T − Tmedel), där Tmedel = 30.7 är medeltemperaturen. Ac = y skrivs nu
1 −10.7 1 −5.2 1 −0.5 1 6.1 1 10.3
c1 c2
=
8.78 8.93 9.06 9.25 9.40
Ytterligare exempel, centrering
Formulera normalekvationerna ATAc = ATy och lös ut c
5 0
0 285.08
c1 c2
=
45.42 8.333
⇒ c =
9.0840 0.0292
Detta ger L0 = c1 − c2 ∗ Tmedel = 8.1866 och
λ = L1/L0 = c2/L0 = 0.0036 vilket överensstämmer med första uträkningen.
Ytterligare exempel, centrering
Avrunda till tre korrekta siffror i normalekvationerna ATAc = ATy
5 0 0 285
c1 c2
=
45.4 8.33
Lös ut c och beräkna λ !
c =
9.0800 0.0292
Detta ger L0 = c1 − c2 ∗ Tmedel = 8.1827 och
λ = L1/L0 = c2/L0 = 0.0036 vilket överensstämmer med första
Matlabkod, ytterligare exempel (centrering)
clear all % Tömmer alla variabler
close all % Stänger alla grafikfönster clc % Tömmer kommandofönstret
x = [20.0 25.5 30.2 36.8 41.0]’; % Temperatur y = [8.78 8.93 9.06 9.25 9.40]’; % Längd
xm = mean(x); % medelvärdet av x
A = [ones(size(x)) (x - xm)];
% operatorn \ löser ut c i minsta kvadraten mening!
c = A\y
residual = A*c - y % beräknar felet i varje mätpunkt
felkvadratsumma = residual’*residual
% alternativt: norm(residual, 2)^2 eller sum(residual.^2)
% Förfina x-intervallet och beräkna kurva
xP = min(x) : (max(x) - min(x))/100 : max(x);
yP = c(1) + c(2)*(xP - xm);
% Plotta mätpunkter och kurva plot(x, y, ’*’, xP, yP)
% Plotta residualen figure
Ickelinjär modell
För att kunna lösa en ickelinjär modell försöker man hitta en linjär ersättningsmodell.
Exempel
Funktionen q = C ∗ pα ska anpassas till en serie mätdata (pi, qi) Den linjära formen får man om man logaritmerar funktionen lnq = lnC + αlnp.
Sätt x = lnp, y = lnq, c1 = lnC och c2 = α.
Problemet har nu omvandlats till ett linjärt problem y = c1 + c2 ∗ x
Efter att det linjära problemet lösts substituerar man tillbaka till ursprungsmodellen.
Residualanalys
Genom att studera hur residualen ser ut kan man få tips om hur funktionen kan korrigeras
Läs i NAM avsnitt 2.5
