1. No category

HH/ITE/BN
Minsta kvadratmetoden och Mathematica
1
Något om Minsta kvadratmetoden
och Mathematica
Bertil Nilsson
2015-08-15
y
y
6
5
xi ,yi
5
Εi
4
4
xi , f xi ,
3
3
2
2
1
0
0
1
2
min
3
4
n
i 1
f xi,
min
5
f x,
6
x
1
yi
yx
2
3
4
5
6
x
2
2
x
ci ,c t mg l
y
6
5
5
4
4
3
3
2
2
1
0
5
10
15
20
25
t s
1
2
3
4
x
2
Minsta kvadratmetoden och Mathematica
HH/ITE/BN
ť Förord
På följande sidor presenteras en elementär "streetwise guide" till minsta kvadratmetoden med flitig användning av Mathematica.
Framställningen är fåordig, fri från pedanteri men i någon mening fullständig. Det man väsentligen behöver veta om begrepp,
terminologi, beteckningar och teori för att modellera och lösa problem i framtida kurser och yrkesliv som ingenjör, naturvetare eller
lärare klarläggs och typiska exempel ges.
ť Inledning
Vi börjar med en problemställning från verkliga livet. Anders på Andersberg har köpt en ny bil med konstantfarthållare. För att
undvika framtida bekymmer med lagens väktare vill nu Anders kontrollera riktigheten i utrustningen. För detta ändamål uppsöker
han en lång raksträcka där han vet att avståndet mellan lyktstolparna är 100 m. Med konstantfarthållaren inställd på 50 km/h bär det
iväg. Eftersom han bara har ett vanligt armbandsur till hands så blir det till att uppskatta restiderna på sekunden när. Han får följande
serie från en stolpe till de fem efterföljande
tid s
0 7 14 22 29 36
Anders strategi är nu att äntligen dra nytta av den gamla mekanikformeln s
börjar med att lägga in mätdata i ett diagram
0 7
14 22 29 36
0 100 200 300 400 500
data
ListPlot data, PlotStyle
vt som han under skoltiden tvingades exercera. Han
;
Red, PointSize 0.03
, AxesLabel
"t
s ", "s
m "
s m
500
400
300
200
100
5
Han förstår att s
10
15
20
25
30
t s
35
vt är en rät linje genom origo, men hur ska den luta? Entusiastiskt låter han den gå genom sista mätpunkten, det
vill säga han väljer modellen s
vt
500
t
36
och lägger till den i diagrammet
500
Plot
t, t, 0, 40 , PlotStyle
Blue, AxesLabel
"t
s ", "s
m " ,
36
Epilog
Red, PointSize 0.03 , Point data

s m
500
400
300
200
100
10
20
30
40
t s
Nja, riktigt bra blev det inte. Genast dyker det upp några frågeställningar. Det kanske finns någon lutning som gör att fler punkter
kommer ännu närmare linjen? Å varför är första mätpunkten 0, 0 helig? Det borde ju vara lika svårt att "pricka" första stolpen som
de andra? I så fall kanske han skulle prova s vt m. Vi ska hjälpa honom att utreda detta.
ť Minsta kvadratmetoden
Antag att vi har givet mätdata och vill anpassa en funktion till dessa "så bra som möjligt". Låt denna funktion allmänt tecknas
y f x, , där x, y är den oberoende respektive beroende variabeln och en vektor med parametrar som ska bestämmas så
anpassningen blir så bra som möjligt. Naturligtvis kan vi ha flera oberoende variabler, t.ex. T f x,
f x, y, z, t ,
som
modell för hur temperaturen varierar med tiden t på olika platser i rummet x, y, z . Vanligtvis kallas f för modellen och för
modellparametrarna. Exempelvis
HH/ITE/BN
Minsta kvadratmetoden och Mathematica
y
kx
y
a
y
m
k, m
bx
c0
3
a, b
2
c1 x
c2 x
c0 , c1 , c2
Det typiska är att vi måste välja modellen själva, se figur nedan där tre olika modeller är anpassade till samma mätdata. Minsta
kvadratmetoden (MKM) hjälper sedan till med att bestämma parametrarna . Vi får alltså olika anpassningar beroende på vilka val
vi gör. Detta är viktigt att förstå! Ofta är det så att man vet vilka funktioner som ska ingå i f. Kanske har man en matematisk modell
där man behöver trimma några parametrar mot försök.
I Minsta kvadratmetoden måste man välja modellen själv och man får olika anpassningar beroende på vilka val man gör!
Det är viktigt att man väljer funktioner i modellen så att den på ett kvalitativt sätt fångar det man vet om systemet som studeras, t.ex.
y
då x 0 eller y 0 då x
. I figuren nedan ser vi också att det är mycket vanskligt att använda modellen utanför mätområdet. Detta kallas extrapolation till skillnad mot interpolation då man befinner sig i det område som omsluts av mätområdet.
Notera också i figuren att med ett komplicerat modellval så kan man lätt få den att gå genom alla mätpunkter, exempelvis ett
polynom av grad n definieras entydigt av n 1 mätpunkter. Men detta kan vara en mycket dålig representation av den bakomliggande fysikaliska modellen, som kanske bättre speglas av de två "mindre nervösa" modellerna, speciellt om man tänkt sig tillåta
"lite" extrapolation. Spridningen i mätdata kan helt enkelt bero på svårigheter vid mätningen. Den kraftiga översvängningen vid
ändpunkterna av mätområdet som uppstår vid approximation med polynom av höga gradtal brukar kallas Runges fenomen, Carl
Runge (1856-1927), se den gröna kurvan.
y
5
4
f x
kx m
3
f x
ax2 bx c
2
f x
c6 x6 c5 x5
c1 x c0
1
1
Interpolation
3
4
2
5
6
x
7
Minsta kvadratmetoden vilar på gamla kända kunskaper, nämligen minimering av en i detta fall kvadratisk funktion. Därav namnet.
Antag att vi har en uppsättning mätdata med n st mätpunkter, se fig nedan från vänster till höger. Bilda nu vid varje mätpunkt xi , yi
och det uppmätta värdet yi , det vill säga Εi f xi ,
yi . Nu är det lika illa om mätpunkten ligger
felet Εi mellan modell f xi ,
över som under modellen. Välj därför i stället att studera Ε2i och lägg samman dem alla.
y
6
xi ,yi
5
x1 y1
x2 y2
Εi
4
xi , f xi ,
S
3
xn yn
n 2
i 1Εi
Ε21
Ε22
Ε2n
n
i 1
f xi ,
yi
2
2
1
0
0
1
2
3
4
5
6
x
En kvadratsumma är alltid icke-negativ, det vill säga S 0, och i det ideala fallet S 0 går modellen genom alla mätpunkter.
Eftersom alla mätpunkter är applicerade blir S en funktion av parametrarna, det vill säga S . Minsta kvadratmetoden går nu ut på
att minimera S med avseende på parametrarna .
min S
Man säger att f är anpassad i minsta kvadratmetodens mening. Vi drar oss nu till minnes att söka minimum till en funktion är det
samma som att söka nollställe till derivatan. Eftersom vi har flera variabler i måste vi simultant söka nollställe till respektive
derivata. Att derivera en funktion av flera variabler med avseende på någon speciell variabel kallas för partiell derivata och skrivs
med i stället för . Inget märkvärdigt, man gör som vanligt, det vill säga deriverar med avseende på den variabel som önskas och
låter de andra vara konstanter vilka som helst. Vi tar ett
Exempel: Låt funktionen g x, y
x3
5y
4xy2
xsin 2 y
7 vara given. Bestäm de partiella derivatornna
g
x
och
g
y
.
4
Minsta kvadratmetoden och Mathematica
Lösningsförslag: Vi får direkt
g
3x2
x
Så att söka minimum av S, min S
4y2
g
sin 2 y och
5
y
8x y
HH/ITE/BN
2xcos 2 y .
S
, blir då att lösa ekvationssystemet
S
p1
0
S
p2
0
Alltså ett ekvationssystem med lika många ekvationer som vi har modellparametrar
p1 , p2 ,
. Om detta blir linjärt har vi
linjär minsta kvadratmetod annars olinjär. En linjär modell har vi alltid om modellen kan skrivas som en skalärprodukt mellan
. Exempelvis har vi modeller som leder till
ingående funktioner och modellparametrarna, det vill säga f x,
f1 , f2 ,
linjär MKM
y
y
kx
kx
y
x ,
x, 1 ,
m
c0
2
c1 x
k
k, m
2
c2 x
1, x, x
,
c0 , c1 , c2
Vid linjär MKM har S endast ett lokalt minimum, som då också är globalt, med entydig lösning som följd. Tryggt! I detta fall finns
till vår hjälp i Mathematica en lättanvänd funktion som heter just "anpassa", Fit[data,f,x]. Den matas med mätdata i form av
.
en matris med två kolonner xi , yi samt en vektor med och levererar sedan f på direkt användbar form
Låt oss nu använda det typiska exemplet y
processen. Först summan av alla kvadratfel
f x,
kx
m, med modellparametrarna
S
n
i 1
kxi
m
yi
k, m , som modell och arbeta igenom
2
Med inre derivata i färskt minne får vi så ekvationssystemet som bestämmer
min S
S
k
S
m
min k,m S
n
i 12
kxi
m
yi xi
n
i 12
kxi
m
yi
0
n
2
i 1 xi
n
i 1 xi
0
n
i 1 xi
n
i 1 xi yi
n
i 1 yi
k
m
n
(1)
Exempel: Studera Anders problem.
Lösningsförslag: Först den modell som går genom origo, det vill säga s
S
6
i 1
72
142
vti
si
S
v
2
6
i 12
vti
vt
si ti
t v
0
v
. Vi får
6
2
i 1 ti
6
i 1 ti si
där
6
2
i 1 ti
6
i 1 ti si
02
0 0
7 100
222
292
14 200
362
22 300
Provkör Fit[data,f,x] på Anders problem. Med s
s1
vt
2866
29 400
t v
36 500
39 700
2866
v
39 700
13.8521
får vi direkt
Fit data, t , t
13.8521 t
Vi känner igen resultatet från handräkningen ovan. Sedan över till hans fundering s
vt
m
t, 1 
v

m
. Vi startar direkt
med slutresultatet av minimeringen (1) ovan. Kontrollräkna gärna!
6
2
i 1 ti
6
i 1 ti
6
i 1 ti
n
v
 
m
6
i 1 ti si
6
i 1 si
2866 108 v
 
m
108
6
39 700
1500
v
 
m
1
461
6350
950
13.7744
2.06074
Visst, Mathematica instämmer
s2
Fit data, t, 1 , t
13.7744 t
2.06074
En liten jämförelse visar att han har mätt väl. Modellerna är inte lika och i detta fall är skillnaden nästan av akademiskt intresse, men
indikerar ändå det viktiga att olika funktionsval ger olika modeller. Vi avslutar med att muppa in data och båda modellerna i samma
diagram. Detta återkommer i tid och otid, bara så du vet !
HH/ITE/BN
Plot
Minsta kvadratmetoden och Mathematica
5
s1 , s2 , t, 0, 40 , PlotStyle
Blue, Orange , AxesLabel
Epilog
Red, PointSize 0.03 , Point data
"t", "si ,s1 ,s2 " ,
si ,s1 ,s2
500
400
300
200
100
10
20
30
40
t
Om vi har en linjär modell kan vi också komma fram till det önskade ekvationssystemet genom att studera det överbestämda
ekvationssystemet som i fallet y f x,
kx m blir
x1 1
x2 1
xn 1
k
m
y1
y2
yn
Detta överbestämda ekvationssystem har i allmänhet ingen lösning eftersom vi har fler ekvationer än obekanta. Man kan visa att det
kan ges en lösning i minsta kvadratmetodens mening, och vi börjar med att övertyga oss genom en geometrisk betraktelse.
är då en vektor i planet med som komponenter i basen k . Eftersom inte löser
Kolonnerna k i spänner upp ett plan och
ligger följaktligen inte i planet. Då är residualvektorn
en vektor som går från spetsen på till spetsen på
i
planet. Minsta kvadratmetoden motsvarar då att söka minimum av
. Vi kommer ihåg att närmsta vägen från en punkt till ett plan
då är vinkelrät mot alla basvektorerna k i planet, det vill säga
går vinkelrät mot planet. Detta innebär att vi har minimum av
kolonnerna i .
Mera formellt har vi min
2
1
0
2
0
min
2
2
.
Vi kallar
för det utjämnade systemet eller normalekvationerna till det överbestämda ekvationssystemet
ovan.
Koefficientmatrisen
är alltid både kvadratisk och symmetrisk beroende på just konstruktionen
. Man kan visa att normalek1
1
, där
kallas pseudoinversen till . Naturligtvis bestäms inte
vationerna alltid har en lösning
på detta kostsamma sätt med inverser utan med Gausselimination direkt på normalekvationerna.
Om
är överbestämt kallas
normalekvationerna och
har entydig lösning i minsta kvadratmetodens mening.
Vi provar på Anders modell s vt m. Vi ser att koefficientmatrisen
är både kvadratisk och symmetrisk som sig bör. Kontrollräkna gärna matrisgodiset och ekvationslösningen som följer, vilken åter igen levererar parametrarna i minsta kvadratmetodens
mening
,1 &
0
7
14
22
29
36
1
1
1
1
1
1
NSolve
m

data All, 1 ;
0 7 14 22 29 36
2866 108
 

1 1 1 1 1 1
108
6
. . v, m
2.06074, v
13.7744
Med samma resultat som ovan.
.data All, 2
.
.data All, 2
39 700, 1500
6
Minsta kvadratmetoden och Mathematica
HH/ITE/BN
Så här via normalekvationerna kan man alltid gå då man har linjär minsta kvadratmetod. Vid olinjär modell är man placerad i osäker
terräng och har man väl hittat ett minimum kan det alltid ifågasättas om det är ett globalt. Till hjälp finns FindFit,
(N)Minimize och (N)Maximize som utger sig för att kunna leverera globalt optimum. De två sistnämnda kan även matas med
bivillkor. Fördelen med FindFit är att man inte behöver formulera kvadratsumman explicit utan kan ange mätdata och den
olinjära modellen som indata. Utöver dessa finns naturligtvis de gamla kända FindMinimum och FindMaximum som då i någon
mening får anses vara något omsprungna av de tre tidigare nämnda. Det finns även specialpaket Statistics`NonlinearFit`, NumericalMath`PolynomialFit`, NumericalMath`SplineFit` och NumericalMath`TrigFit`.
Hittills har vi behandlat diskret minsta kvadratmetod eftersom vi anpassat en modell till ett ändligt antal mätpunkter. Om vi i någon
mening låter antalet mätpunkter växa obegränsat får vi kontinuerlig minsta kvadratmetod.
min S
min
n
i 1
f xi ,
2
yi
min
f x,
yx
2
x
n
Som sinnebild kan man se detta som att man minimerar den målade arean som innestängs mellan (modell)funktionen f x,
och
funktionen y x i intervallet , fast inte riktigt eftersom vi kvadrerar integranden. Detta kan komma till användning då man i ett
intervall vill approximera en besvärlig funktion med någon enklare. Även här har vi de två fallen med linjär och olinjär minsta
kvadratmetod.
Exempel: Som illustration till kontinuerlig minsta kvadratmetod väljer vi att i intervallet
2
med den betydligt enklare linjära modellen f x,
och
a, b, c de sökta modellparametrarna.
ax
bx
2
c
x , x, 1
a, b, c
0, 1 approximera y
, där
x
1
sin
10
10x
2
x , x, 1 är funktionsvalen
Lösningsförslag: Med ett snabbt hack och ointressanta utdataceller undertryckta
1
a, b, c ; f
x2 , x , 1. ; S

f
0
x
2
1
Sin 10 x
x;
10
får vi så en inblick i skådespelet
Plot
x
1
Sin 10 x , f . Solve D S,
0&
First, x, 0, 1 ,
10
PlotStyle
Red, Blue , AxesLabel
"x", "y x ,f x,
" 
y x , f x,
2.5
2.0
1.5
0.2
0.4
0.6
0.8
1.0
x
Som läsaren säkert upptäckt nu så är minsta kvadratmetoden inget som räknas för hand. Eftersom den är mycket vanlig i praktiken
och då kanske matas med flera tusen mätpunkter och modeller med många parametrar och komplicerade funktioner överlämnas
gärna råräknandet till datorerna. Det viktiga är att man vet vad man gör innan man lämnar över till Mathematica som har effektiva
lösare för alla förekommande minsta kvadratmodeller. Som avslutning ska man rita modell och mätdata i samma diagram och göra
en visuell betraktelse av hur noggrant modellen ansluter till mätdata. Detta är inte helt lätt då antalet oberoende variabler i modellen
är större än två, då får man studera olika modeller och deras kvadratfel. Med andra ord så ligger det lite "godtycke" med i metoden
och vid komplicerade problem bör någon erfaren person konsulteras så man inte sprider illäror med en modell som är helt uppåt
väggarna!
Var noga med modellval, "man får den anpassning man förtjänar" och undvik interpolation i mätområdets utkanter
samt befatta dig inte med extrapolation!
Resten av skriften är en bukett exempel, de flesta utan bakgrund och med få mätpunkter för att hålla beräkningsvolymen på en
beskedlig nivå. Utanför skolans värld krävs definitivt information om vilket system som det är mätt på och fler mätpunkter för att
kunna göra ett välunderbyggt val av modell! Vi avslutar med några sådana exempel.
HH/ITE/BN
Minsta kvadratmetoden och Mathematica
7
ť Blandade exempel för resten
Exempel: Använd MKM för att anpassa en rät linje p1 x
spolynom p3 x c3 x3 c2 x2 c1 x c0 till mätpunkterna
x
kx
m, ett andragradspolynom p2 x
1
2
ax2
bx
c och ett tredjegrad-
3.5 5.8
y 1.1 1.8
4
5.5
Lösningsförslag: Typisk råräkning. Börja med räta linjen och möblera det överbestämda ekvationssystemet
x1 1
x2 1
k
m
xn 1
y1
y2
yn
med hjälp av givna mätdata
1
2 3.5 5.8
1.1 1.8 4 5.5
data
,1 &
1
2
3.5
5.8
;
data All, 1
1
1
1
1
Bestäm nu k och m enligt MKM receptet och normalekvationerna
.
. Vi tar lite matrisgodis på vägen
.data All, 2
1 2 3.5 5.8
50.89 12.3
 

1 1 1 1
12.3
4.

k
m
50.6, 12.4
Slutligen de efterlängtade.
Solve
k
. . k, m
0.954276, m
.data All, 2
0.165602
Å så en liten bild
Plot k x m .
Epilog
, x, 0, 7 , PlotStyle Blue, AxesLabel
Red, PointSize 0.03 , Point data
"x", "yi ,f x,
" ,
yi , f x,
7
6
5
4
3
2
1
1
2
3
4
5
6
7
x
Det går även bra att direkt använda funktionen Fit. Här väljer man vilka funktioner man vill använda, i retur kommer en linjärkombination av dessa i MKM:s mening.
p1
Fit data, 1, x , x
0.954276 x
0.165602
Med samma resultat som ovan. Vi provar nu med ett andragradspolynom, p2 x
att möblera det överbestämda ekvationssystemet
ax2
bx
c. Detta går naturligtvis att lösa genom
8
Minsta kvadratmetoden och Mathematica
x21 x1 1
x22 x2 1
x2n xn 1
Först
a
b
c
y1
y2
yn
, sedan lite matrisgodis och slutligen parametrarna ur normalekvationerna
2

, , 1 &
HH/ITE/BN
.
data All, 1
1
1 1
4
2 1
12.25 3.5 1
33.64 5.8 1
.
.data All, 2
1 4 12.25 33.64
1 2 3.5
5.8
1 1
1
1
Solve
a
1298.71 246.987 50.89
246.987 50.89
12.3
50.89
12.3
4.
. . a, b, c
0.0735314, b
242.32, 50.6, 12.4
.data All, 2
1.46352, c
0.464835
Man ska kunna genomföra kalkylen den hårda vägen via normalekvationerna ovan, men att för hand multiplicera ihop och lösa ett
ekvationssytem med tre obekanta är direkt olämpligt i datorernas tidevarv! Det viktiga är att förstå hur det fungerar, sedan överlämnar man gärna råräknandet! I praktiken hade man valt den lätta vägen med Fit[data,f,x]och istället lagt kraft på att prova
några olika modeller
Fitdata, 1, x, x2 , x
p2
0.0735314 x2
1.46352 x
0.464835
Avslutningsvis de två modellerna tillsammans med mätdata. Jag har varnat
Plot
p1 , p2 , x, 0, 7 , PlotStyle
Blue, Orange ,
AxesLabel
"x", "yi ,p1 ,p2 " , Epilog
Red, PointSize 0.03 , Point data
yi , p1 , p2
6
4
2
1
2
3
4
5
6
7
x
Avslutningsvis modellen med tredjegradspolynom direkt med Fit[data,f,x].
Fitdata, 1, x, x2 , x3 , x
p3
0.108543 x3
1.0122 x2
1.57679 x
1.77314
Det är nu bara att välja den modell som önskas beware
Plot
p1 , p2 , p3 , x, 0, 7 , PlotStyle
Blue, Orange, Green ,
AxesLabel
"x", "yi ,p1 ,p2 ,p3 " , Epilog
Red, PointSize 0.03 , Point data
yi , p1 , p2 , p3
6
4
2
1
2
3
4
5
6
7
x
HH/ITE/BN
Minsta kvadratmetoden och Mathematica
Exempel: Anpassa med minsta kvadratmetoden y x
9
bx2 till mätvärdena
a
x 5
7.5 12 15 25
y 13.1 18.1 70.2 109 301
Lösningsförslag: De fem mätvärdena möblerar
och
i det överbestämda ekvationssystemet
1 x21
a
b
x22
1
y1
y2
a
b
för de sökta parametrarna a och b, där
5
7.5
12
15 25
13.1 18.1 70.2 109 301
data
2
1,
&
;
data All, 1
1 25
1 56.25
1 144
1 225
1 625
a
b
De efterlängtade parametrarna a och b levereras nu av normalekvationerna
.

.data All, 2
1
1
1
1
1
5.
1075.25
 

25 56.25 144 225 625
1075.25 465 775.
Solve
a
. Kontrollräkna!
. . a, b
2.36283, b
511.4, 224 104.
.data All, 2
0.486598
Eller direkt med Fit.
Fitdata, 1, x2 , x
f
0.486598 x2
2.36283
En bild piggar alltid upp
Plot f, x, 0, 25 , PlotStyle Blue, AxesLabel
Epilog
Red, PointSize 0.03 , Point data
"x", "yi ,f x,
" ,
yi , f x,
300
250
200
150
100
50
5
10
15
20
25
x
Exempel: Anpassa med minsta kvadratmetoden y x
ax
bx2 till mätvärdena
x 2
3
y 1.9 0.1
Lösningsförslag: De fyra mätvärdena möblerar
och
x1 x21
x2
x22
5
11
8
39
i det överbestämda ekvationssystemet
a
b
y1
y2
a
b
10
Minsta kvadratmetoden och Mathematica
HH/ITE/BN
för de sökta parametrarna a och b, där
2
3
1.9 0.1
data
2
 ,
&
5
11
8
39
;
data All, 1
2 4
3 9
5 25
8 64
Bestäm nu a och b enligt MKM receptet och normalekvationerna
.

. Kontrollräkna!
.data All, 2
2 3 5 8
102 672
 

4 9 25 64
672 4818
Solve
a
a
b
. . a, b
2.70872, b
362.9, 2762.5
.data All, 2
0.951174
Eller direkt med Fit.
Fitdata, x, x2 , x
f
0.951174 x2
2.70872 x
Javisst, vi vill se en bild
Plot f, x, 2, 8 , PlotStyle Blue, AxesLabel
Epilog
Red, PointSize 0.03 , Point data
"x", "yi ,f x,
" ,
yi , f x,
3
4
5
6
7
8
x
10
20
30
40
Exempel: Kväve och syre har atommassorna N 14 och O 16. I tabellen nedan finns experimentellt uppmätta molekylmassor för
sex olika kväveoxider. Använd minsta kvadratmetoden för att bestämma kvävets och syrets atommassor.
NO
N2 O NO2 N2 O3 N2 O4 N2 O5
30.01 44.01 46.03 75.98 92.02 108.01
Lösningsförslag: Lite olika blandningar av "äpplen och päron" på vågen. De sex vägningarna möblerar
mN

för de sökta atommassorna.
ekvationssystemet 
mO
1 2 1 2 2 2
1 1 2 3 4 5
;
30.01, 44.01, 46.03, 75.98, 92.02, 108.01 ;
Bestäm nu mN och mO enligt MKM receptet och normalekvationerna
.

.
1 2 1 2 2 2
18 29
 

1 1 2 3 4 5
29 56
Solve
mN
. . m N , mO
14.0008, mO
16.0023
716.08, 1302.15
.

mN

mO
. Kontrollräkna!
och
i det överbestämda
HH/ITE/BN
Minsta kvadratmetoden och Mathematica
Exempel: Anpassa med minsta kvadratmetoden y x
x
ax b
11
till mätvärdena
x 1 2.5 4 6 8
y 1 1.5 2 2.5 3
Lösningsförslag: Här har vi en olinjär modell, men kan efter omskrivning lineariseras och möblera ett överbestämt ekvationssystem.
x
1
ax b x 0 y
y
ax b
x
1
y
a
b
x
1
1
x1
1
1
x2
a
b
1
y1
a
b
1
y2
för de sökta parametrarna a och b, där
1 2.5 4 6 8
1 1.5 2 2.5 3
data
;
1
&
1,
data All, 1
1 1
1 0.4
1
1
4
1
1
6
1
1
8
Bestäm nu a och b enligt MKM. Kontrollräkna gärna matrisgodiset och ekvationslösningen som följer
.

1
1
1
.
data All,2
1 1 1
1 0.4
1
1
1
4
6
8


5.
1.94167

1.94167 1.2659
2.9, 1.5
1
Solve
. . a, b
.

data All, 2
a
0.2964, b
0.730302
Här kan det vara lämpligt att introducera den olinjära hjälpredan FindFit som sväljer modellen med hull och hår
x
, a, b , x
FindFitdata,
ax
a
0.212872, b
b
1.0605
Ok, dé får väl bli en bild då
x
Plot
.
ax
, x, 1, 9 , PlotStyle
PlotRange
0, 4 , Epilog
"x", "yi ,f x,
Red, PointSize 0.03 , Point data
yi , f x,
4
3
2
1
2
Blue, AxesLabel
b
4
6
8
x

" ,
12
Minsta kvadratmetoden och Mathematica
Exempel: Sambandet mellan trycket p och volymen v i en gas följer lagen pvn
minsta kvadratmetoden och mätserien
HH/ITE/BN
c, där n och c är konstanter. Bestäm dessa med
v 4.60 7.20 10.1 15.3 20.4 30.0
p 14.2 7.59 4.74 2.66 1.78 1.04
Lösningsförslag: Här har vi en olinjär modell, men kan lineariseras efter logaritmering av båda sidor
n
pv
c
ln p
nln v
ln c
ln c
nln v
ln p
1
1
ln v1
ln v2
ln c
n
ln p1
ln p2
ln c
n
Nu finns det fotspår i sanden
4.60 7.20 10.1 15.3 20.4 30.0
14.2 7.59 4.74 2.66 1.78 1.04
data
1,
1
1
1
1
1
1
Log
&
;
data All, 1
1.52606
1.97408
2.31254
2.72785
3.01553
3.4012
Bestäm nu ln c och n enligt MKM receptet och normalekvationerna
.

lnc
.
.Log data All, 2
1
1.52606
lncÅn
ln c
n
1
1.97408
Solve
4.77908, n
1
2.31254
1
2.72785
. . lnc, n
1
3.01553
1
6.
 
3.4012
14.9573
.Log data All, 2
14.9573

39.6764
First
1.39359
Så modellen
lnc
p
. lncÅn
vn
118.995
v1.39359
Eller direkt med den olinjära hjälpredan
c
FindFitdata,
, c, n , v
vn
c
119.337, n
1.39505
Detta eviga ritande
Plot p, v, 0, 25 , PlotStyle Blue, AxesLabel
Epilog
Red, PointSize 0.03 , Point data
pi , p v
35
30
25
20
15
10
5
5
10
15
20
25
v
"v", "pi ,p v " ,
7.83027, 16.1894
HH/ITE/BN
Minsta kvadratmetoden och Mathematica
Exempel: Efter att en sjö försetts med ett antal fiskyngel uppskattades,
bl.a. från åtgången av mat, hur antalet fiskar f utvecklades över tiden t
mätt i år. Man fick följande tabell
t
f
1
2
3
4
5
6
7
8
9
10
fiskar
1000
800
.
100 210 360 550 720 840 910 940 950 960
600
400
I många biologiska sammanhang visar det sig att den så kallade logistiska
f
modellen fungerar mycket bra. Denna har formen f t
1 
f
f0
,
kt
1
13
200
där de tre parametrarna f , f0 och k styr dess utseende. Anpassa denna
0
2
4
6
8
10
t år
modell till data i tabellen och uppskatta hur många fiskyngel som planterades ut från början f0 samt hur många fiskar det kommer
att finnas vid jämvikt f , det vill säga efter lång tid.
Lösningsförslag: Vi har en typisk ickelinjär modell, med data och föreslagen modell.
f
1
2
3
4
5
6
7
8
9
10
100 210 360 550 720 840 910 940 950 960
data
; modell
;
1

f
kt
1
f0
Vi tar direkt hjälp av Mathematicas funktion för olinjär MKM. Vi kan därefter direkt avläsa de önskade parametrarna och rita en
liten figur ända från år noll då utplanteringen skedde.
FindFit data, modell, f , f0 , k , t
967.146, f0
f
49.6372, k
0.799629
Plot modell . , t, 0, 10 , AxesLabel
"t
PlotRange
Automatic, 0, 1000 , Epilog
år ", "fiskar" , PlotStyle Blue,
Red, PointSize 0.03 , Point data
fiskar
1000
800
600
400
200
0
2
4
6
8
10
t år
Exempel: Inte sällan har man nödvändigtvis en modell som genererar ett icke linjärt system. Ta som exempel tomtegröt som ju
avsvalnar enligt Newtons avsvalningslag
T
t
k T0
T . Lösningen får den välkända formen T t
a
kt
b
, där parametrarna a, b
och k ingår olinjärt. Antag att tomten nedtecknat följande mätserie.
t min
1
3
5
10 15 20
T C
75 55 40 28 25 23
n 2
i 1Εi .
Lösningsförslag: På grund av olinjäriteten måste vi angripa problemet med en direkt minimering av felsumman S
Mätdata
och modell
1 3 5 10 15 20
75 55 40 28 25 23
data
;T t
: a
kt
b
Felet vid varje mätpunkt
Ε
 a
data All, 2
b
Felsumman S
S
 a
k
75, a
n
2
i 1 Εi
T data All, 1
b
3k
55, a
b
5k
40, a
b
10 k
28, a
b
15 k
25, a
b
20 k
23
får vi nu enklast om vi betraktar Ε som en vektor. En enkel skalärprodukt gör jobbet.
Ε.Ε
b
20 k
2
23
 a
b
15 k
2
25
 a
b
10 k
2
28
 a
b
5k
2
40
 a
Sedan minimering. FindMinimum kräver som vanligt en startpunkt där den ska börja leta
b
3k
2
55
 a
b
k
2
75
14
Minsta kvadratmetoden och Mathematica
Smin ,
HH/ITE/BN
FindMinimum S, a, 20 , b, 50 , k, 0.5
3.4631, a
23.0881, b
68.0623, k
0.264998
Alltså modellen
T t
.
23.0881
68.0623
0.264998 t
Visst blir man glad av upprepningar
Plot T t
. , t, 0, 30 , PlotStyle Blue, AxesLabel
"t min ", "Ti ,T t
PlotRange
0, 100 , Epilog
Red, PointSize 0.03 , Point data
Ti ,T t
100
C " ,
C
80
60
40
20
0
5
10
15
20
25
30
t min
Om vi namnmässigt korsar FindMinimum med Fit får vi den i särklass enklaste funktionen att använda. Inte ens startvärde
behövs i det olinjära fallet!
FindFit data, T t , a, b, k , t
a
23.0881, b
68.0623, k
0.264998
Exempel: För att studera en människas rörelser under gång har man utvecklat en datormodell. För att generera indata till denna
utgår man från fotografier som tagits vid åtta tillfällen under två steg, en period. Tiden för en period är normerad till 1 för att enklare
kunna simulera olika gånghastigheter. I nedanstående tabell är knävinkeln i grader angivet för de olika tidpunkterna. För att unvika
minst sagt ryckig gång på datorskärmen vill man givetvis anpassa en funktion som levererar en vinkel vid godtycklig tidpunkt t.
data
0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1
180 160
165
175 170 130
110
150 180
;
Lösningsförslag: Här är det viktigt att använda en periodisk funktion så att vi får en kontinuerlig gånghastighet (derivata) över
periodgränsen. Använd en harmonisk funktion som automatiskt tar hand om periodiciteten, en så kallad Fourierserie, här
3
t
bk sin 2Πkt , med modellparametrarna
a0 , a1 , b1 , a2 , b2 , a3 , b3 eftersom b0 0. Så modellen direkt
k 0 ak cos 2Πkt
färdigmöblerad.
f
1, Cos 2 Π
t , Sin 2 Π
t
&
1, 2, 3
Flatten
1, cos 2 Π t , sin 2 Π t , cos 4 Π t , sin 4 Π t , cos 6 Π t , sin 6 Π t
Fit data,
23.4727 sin 2 Π t
Varför inte rita lite
f,
t
8.75 sin 4 Π t
4.02728 sin 6 Π t
3.55055 cos 2 Π t
18.9167 cos 4 Π t
1.78278 cos 6 Π t
155.083
HH/ITE/BN
Plot
Minsta kvadratmetoden och Mathematica
15
, t, 0, 1 , PlotStyle Blue, AxesLabel
"t", " i , t, " ,
PlotRange
100, 200 , Epilog
Red, PointSize 0.03 , Point data
i , t,
200
180
160
140
120
0.0
0.2
0.4
0.6
0.8
1.0
t
Eftersom det i båda ändpunkterna är uppmätt 180 bör man kan kanske kontrollera att ingen översvängning skett. Knävinklar över
180 är klart ohälsosamma
FindMaximum
179.523, t
, t, 1
0.99025
Exempel: En teknik som används på sjukhus för att mäta kapaciteten (volym och blodflöde) i ett inre organ i kroppen är att använda
sig av ett spårämne eller kontrastvätska. En känd mängd av spårämnet injiceras i ingående artär strax intill organet. I utgående ven
mäter man sedan hur spårämnets koncentration c i blodet varierar med tiden. I tabellen nedan ses resultatet av en sådan undersökning sedan 5 mg av ett spårämne injicerats.
tid s
0
1
2
3
4
c mg l
0
0
0
0.1
1.8
tid s
5
6
7
8
9
c mg l
3.0
3.8
4.4
4.9
5.2
tid s
10
11
12
13
14
c mg l
5.4
5.6
5.7
5.8
4.2
tid s
15
16
17
18
19
c mg l
3.0
2.1
1.5
1.1
0.8
tid s
20
21
22
23
24
c mg l
0.6
0.4
0.4
0.2
0.1
Uppgiften är nu att ta fram en modell där såväl blodflödet i [l/s] som organets volym i [l] kan bestämmas.
Lösningsförslag: Mätdata enligt uppgift. Tid i kolonn ett och koncentration i kolonn två.
data
0 0 5
1 0 6
SortPartitionFlatten 2 0 7
3 0.1 8
4 1.8 9
3.0
3.8
4.4
4.9
5.2
10
11
12
13
14
5.4
5.6
5.7
5.8
4.2
15
16
17
18
19
3.0
2.1
1.5
1.1
0.8
20
21
22
23
24
0.6
0.4
0.4 , 2;
0.2
0.1
En bild förtydligar det som berättas i tabellen.
ListPlot data, Joined True, AxesLabel
"t s ", "ci mg l " ,
PlotStyle Red, Epilog
Red, PointSize 0.025 , Point data
ci mg l
6
5
4
3
2
1
5
10
15
20
t s
När koncentrationen ökar från noll till dess den klingar av från maximum, det vill säga i intervallet [3,13] s, har vi själva injiceringsfasen på andra sidan organet, alltså under 10 s. Antag att vi har ett organ med konstant volym och perfekt omrörning, samt att
injektionen sker med konstant flöde under de 10 sekunderna. Från en kurs i analys känner vi kanske igen differentialekvationen för
kontinuerlig blandning
16
Minsta kvadratmetoden och Mathematica
q,cin
q,ct
V , ct
HH/ITE/BN
c
t
V
q cin
c
där c c t [mg/l] är koncentrationen av spårämnet i organet samt de sökta modellparametrarna V [l] organets volym, q [l/s] flödet
genom organet och cin [mg/l] inkommande spårämnets koncentration. Tillsammans med begynnelsevärdet c 3 0.1 blir lösningen
c t till begynnelsevärdesproblemet
1
cAvt
DSolveV c ' t
q Cin
c t
,c 3
, c t , t
First
FullSimplify
10
1
q t 3
c t
1
10 Cin
10
Cin 
V
Denna icke-linjära modell ska nu anpassas till mätdata, under själva injiceringsfasen 3, 13 s, med hjälp av minsta kvadratmetoden.
T, CT

data Range 4, 14
3
4 5 6 7
8
9 10 11 12 13

0.1 1.8 3. 3.8 4.4 4.9 5.2 5.4 5.6 5.7 5.8
Felet i varje mätpunkt är skillnaden mellan det uppmätta värdet och vad modellen ger.
CT
Ε
c t
. cAvt . t
1
 Cin
T
1
10
10 Cin
1
q
0.1,
1
10
1
10
1
10 Cin
V
Cin
3.,
10
1
10 Cin
V
Cin
4.9,
10
10 Cin
V
Cin
5.6,
1
Cin
V
10 Cin
10
1
V
1
10 Cin
10
V
10 Cin
10
1
Cin
5.2,
Cin
5.7,
Cin
V
4.4,
7q
1
10 Cin
1
10 Cin
10
1
9q
1
4q
3.8,
6q
1
2
i Εi .
Bilda kvadratsumman av felen, S
10 Cin
10
1
8q
1
1.8,
3q
1
5q
1
Cin
V
1
2q
1
10 Cin
V
Cin
5.4,
Cin
5.8
10 q
10
V
Beräknas lättast med skalärprodukt. Lägg även till kravet att hela dosen om 5 mg ska vara
på plats efter injiceringsfasen 3, 13 s, det vill säga under 10 s.
S
Ε.Ε
Cin q 10
1
5
2
2
10 q
1
10 Cin
10
Cin
V
1
10
10 Cin
V
1
10
10
10 Cin
V
1
Cin
10
10 Cin
V
Cin
10
10 Cin
Cin
V
2
Cin
1.8
10 q Cin
5
2
1
5.2
10
1
10 Cin
5.6
V
2
Cin
4.9
Cin
3.
2
2q
10
V
2
1
Cin
10 Cin
1
3.8
Cin
V
5q
1
10
3q
1
10 Cin
10
2
2
8q
1
2
V
1
4.4
q
1
10 Cin
1
5.7
6q
1
2
Cin
V
1
5.4
4q
1
10 Cin
10
Cin
2
9q
1
2
7q
1
1
5.8
10 Cin
1
0.1
Vi har ett icke-linjärt minsta kvadratproblem, det vill säga de sökta parametrarna V , q och cin kan inte bestämmas ur ett linjärt
ekvationssystem. På grund av doseringskravet om 5 mg kan vi inte använda den olinjära hjälpredan FindFit, så vi får tillgripa en
generell minimeringsalgoritm
Smin ,
FindMinimum S, V, 0.1 , q, 0.1 , Cin , 6
0.00555799, V
0.250346, q
0.0833842, Cin
5.99634
Varav den anpassade modellen, som naturligtvis bara gäller för t
modell
c t
cAvt .
5.99634
5.89634
0.333076 t 3

En enkel jämförelse fås genom att lägga kurvorna ovanpå varandra.
3, 13 s.
HH/ITE/BN
Minsta kvadratmetoden och Mathematica
17
Plot c t
. modell, t, 3, 13 , AxesLabel
"t s ", "ci ,c t
mg l " ,
PlotRange
0, 25 , Automatic , PlotStyle
Blue, Thickness 0.02 ,
Epilog
Red, Line data , PointSize 0.03 , Point data
ci ,c t mg l
6
5
4
3
2
1
0
5
10
15
20
25
t s
Ser bra ut. God överensstämmelse i det intervall, 3, 13 s, som anpassning skett. Även kravet på den injicerade mängden om 5 mg
är uppfyllt
Cin q 10
5
7.88861 10
31
2
.
Mot bakgrund av ovanstående är det rimligt att svara:
Organets volym V 250 ml, blodflöde q 83 ml/s och koncentration under injiceringsfasen cin
6 mg/l utom tävlan.
Exempel: I många industriella tillämpningar behöver man mäta upp cirklar eller
cirkelbågar. Ett utmärkt exempel är SKF med sina kul och rullager. Utifrån en
given mätserie med x, y –punkter vill man få reda på centrum och radie på den
cirkel som bäst anpassar sig, samt vilken noggrannhet man har. Här kommer
minsta kvadratmetoden till användning.
Lösningsförslag: En mätserie med x, y -punkter från en uppmätning är givna. En bild förtydligar som vanligt situationen.
Graphics
Red, Line Append data, First data
, AxesLabel
"x", "y" , Axes
True
y
5
4
3
2
1
2
3
x
4
Som modell ska vi anpassa en cirkel x
centrum samt dess radie.
x0
2
y
y0
2
r2 , där modellparametrarna x0 , y0 och r är koordinaterna för cirkelns
Vi tycker oss ha ett icke-linjärt minsta kvadratproblem, det vill säga de obekanta modellparametrarna kan inte bestämmas ur ett
linjärt ekvationssystem, så vi tillgriper en generell minimeringsalgoritm för att minimera S. Bilda först felet i varje mätpunkt
Ε
x
2
0
y
0
2
r2
. x
1 ,y
2
&
data;
Eftersom de är ganska många visar vi bara de två första.
Take Ε, 2
 4.07418
0
2
3.1305
Bilda kvadratsumman av felen, S
S
0
2
r2 , 4.06494
2
i Εi .
0
2
3.26086
0
2
r2 
Beräknas lättast med skalärprodukt.
Ε.Ε;
Sök slutligen minimum. Bra startvärde får man genom att snegla på bilden ovan eller en enkel beräkning av medelvärden.
18
Minsta kvadratmetoden och Mathematica
Smin ,
FindMinimum S,
1.77549,
1.99736,
0
0,
2.99473, r
0
2 ,
0,
HH/ITE/BN
3 , r, 2
2.0574
Nu kan vi inspektera skådespelet
Show Graphics
Red, Line Append data, First data
Blue, Circle
.
,
0,
0 , r
AspectRatio Automatic, AxesLabel
,
"x", "y" , Axes
True
y
5
4
3
2
1
2
3
x
4
Ser ju inte så illa ut! Felet i varje mätpunkt, eller åtminstone de fem första
Take fel
Ε .
,5
0.0987187, 0.112806, 0.221534, 0.111521, 0.019211
Till beloppet största felet
Max Abs fel
0.272582
Ofta använt mått är det kvadratiska medelfelet
fel.fel
Length data
0.133247
n 2
i 1Εi
En alternativ formulering får vi genom att som ovan låta S
med Εi
xi
x0
2
yi
y0
2
r2 . Efter derivation med
avseende på x0 och y0 får vi de linjära normalekvationerna (visa gärna detta ;-)
1
n
n
2
i 1 xi
1
n
n
i 1 xi yi
n
2
i 1 xi
n
i 1 xi
1
n
n
i 1 xi yi
n
i 1 yi
n
2
i 1 yi
n
i 1 xi
1
n
n
i 1 yi

n
2
i 1 yi
x0

y0
1
2
n
2
i 1 xi xi
y2i 
n
2
i 1 yi xi
y2i 
1
n
1
n
n
i 1 xi

n
2
i 1xi
y2i 
n
i 1 yi

n
2
i 1xi
y2i 
och slutligen radien
1
n
r
n
i 1
xi
x0
2
yi
y0
2
En tredje variant som också leder till en linjär modell är att utveckla vänsterledet i cirkelns ekvation
x
Låt u
r2
x20
x0
2
y
y0
2
r2
x2
2xx0
x20
y2
2 y y0
y20
r2
2xx0
2 y y0
r2
y20 så har vi efter applicering av mätpunkterna ett linjärt överbestämt ekvationssysten
2x1 2 y1 1
2x2 2 y2 1
Därmed är vi trygga, så låt oss göra kalkylen.
x0
y0
u
x21
y21
x22
y22
x0
y0
u
x20
y20
x2
y2
HH/ITE/BN
Minsta kvadratmetoden och Mathematica
Append 2 , 1 &
. &
data;
19
data;
Sedan lösning av normalekvationerna.
Solve
u
. . x0 , y0 , u
1.99736, y0
8.72498, x0
.
First
2.99473
Varav till slut radien
x20
u
y20
.
2.0574
Avslutningsvis låter vi Mathematica komma till tals med den riktigt snabba varianten ;-)
FindFitAppend
0
2.00759,
0
,0 &
3.00478, r
data, x
2.06441
0
2
y
0
2
r2 ,
0,
0,
r , x, y 