Institutionen för systemteknik
Department of Electrical Engineering
Examensarbete
Adaptiva metoder för systemidentifiering med
inriktning mot direkt viktoptimering
Examensarbete utfört i Reglerteknik vid Tekniska högskolan i Linköping
av
Tony Gillberg
LiTH-ISY-EX--10/4019--SE
Linköping 2010
Department of Electrical Engineering Linköpings tekniska högskola
Linköpings universitet Linköpings universitet
inriktning mot direkt viktoptimering
Examensarbete utfört i Reglerteknik
vid Tekniska högskolan i Linköping
av
Tony Gillberg
LiTH-ISY-EX--10/4019--SE
Handledare: Henrik Ohlsson
isy, Linköpings universitet
Examinator: Martin Enqvist
isy, Linköpings universitet
Division of Automatic Control Department of Electrical Engineering Linköpings universitet
SE-581 83 Linköping, Sweden
2010-06-14 Språk Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport
URL för elektronisk version
http://www.control.isy.liu.se http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-4019 ISBN — ISRN LiTH-ISY-EX--10/4019--SE Serietitel och serienummer Title of series, numbering
ISSN —
Titel Title
Adaptiva metoder för systemidentifiering med inriktning mot direkt viktoptime-ring
Adaptive Bandwidth Selection for Nonlinear System Identification with Focus on Direct Weight Optimization
Författare Author
Tony Gillberg
Sammanfattning Abstract
Direct Weight Optimization (DWO) is a nonparametric system identification method. In DWO the value of a function in a certain point is estimated by a weighted sum of measured values. The weights are obtained as a solution to a convex optimization problem.
DWO has a design parameter which has to be chosen or estimated a priori. There are many ways to estimate this parameter. The main focus of this thesis is to estimate this parameter locally. The advantage of estimating the parameter locally is that the estimate will adapt if the system changes behavior from slowly varying to rapidly varying. Estimation methods of this type are usually called adaptive estimation methods.
There are a number of adaptive estimation methods and the adaptation of some of these methods to DWO has already been done. There are however no evaluation studies done. The goal with this thesis is therefore to find out how to estimate the parameter in DWO in the best way and to find out whether DWO is a good base for an adaptive method.
It turned out that DWO might be too sensitive to local changes in the design parameter to be a good base for an adaptive method. However, one of the adaptive estimation methods stands out from the rest because it is much better than the other methods when it, perhaps, should not. Why this method is good might be a good subject for further research.
Nyckelord
Keywords Direct Weight Optimization, Adaptive Bandwidth Selection, Nonlinear System Identification, Nonparametric System Identification, Local Parameter Estimation
Abstract
Direct Weight Optimization (DWO) is a nonparametric system identification meth-od. In DWO the value of a function in a certain point is estimated by a weighted sum of measured values. The weights are obtained as a solution to a convex optimization problem.
DWO has a design parameter which has to be chosen or estimated a priori. There are many ways to estimate this parameter. The main focus of this thesis is to estimate this parameter locally. The advantage of estimating the parameter locally is that the estimate will adapt if the system changes behavior from slowly varying to rapidly varying. Estimation methods of this type are usually called adaptive estimation methods.
There are a number of adaptive estimation methods and the adaptation of some of these methods to DWO has already been done. There are however no evaluation studies done. The goal with this thesis is therefore to find out how to estimate the parameter in DWO in the best way and to find out whether DWO is a good base for an adaptive method.
It turned out that DWO might be too sensitive to local changes in the design parameter to be a good base for an adaptive method. However, one of the adaptive estimation methods stands out from the rest because it is much better than the other methods when it, perhaps, should not. Why this method is good might be a good subject for further research.
Sammanfattning
Direkt viktoptimering (Direct Weight Optimization, DWO) är en ickeparamterisk systemidentifieringsmetod. DWO bygger på att man skattar ett funktionsvärde i en viss punkt genom en viktad summa av mätvärden, där vikterna optimeras fram. Det faktum att DWO har en inparameter som man måste veta i förväg leder till att man på något sätt vill skatta denna inparameter. Det finns många sätt man kan göra denna skattning på men det centrala i denna uppsats är att skatta inparametern lokalt. Fördelen med detta är att metoden anpassar sig om till exempel systemet ändrar beteende från att variera långsamt till att variera snabbare. Denna typ av metoder brukar kallas adaptiva metoder.
Det finns flera metoder för att skatta en inparameter lokalt och anpassningen till DWO är redan klar för ett fåtal som lämpar sig bra. Det är dock inte undersökt vilken av dessa metoder som ger det bästa resultatet för just DWO. Syftet med denna uppsats är alltså att ta reda på hur man lokalt kan skatta en inparameter
till DWO på bästa sätt och om DWO är en bra grund att basera en adaptiv metod på.
Det har visat sig att DWO kanske är för känslig för en lokalt vald inparameter för att vara en bra grund att basera en adaptiv metod på. Däremot utmärker sig en av metoderna för att skatta inparametern genom att vara mycket bättre än de andra metoderna när den kanske inte borde vara det. Varför den är så bra kan vara ett bra ämne för vidare forskning.
Tack
Först och främst vill jag tacka min handledare Henrik Ohlsson för det oerhörda tålamod som måste ha krävts för att vara med på hela den långa resan och allt annat han hjälpt till med. Tack också till min examinator Martin Enqvist för att han tog på sig ansvaret att examinera detta examensarbete och för värdefulla synpunkter och idéer. Ett speciellt tack till Jacob Roll för idén till examensarbetet och för att han varit tänkt examinator fram till hans jobbyte.
Till sist vill jag tacka min familj och mina vänner och bekanta för allt stöd och all förståelse som har hjälpt mycket.
Innehåll
1 Inledning 1 1.1 Bakgrund . . . 1 1.2 Syfte . . . 2 1.3 Disposition . . . 2 2 Systemidentifiering 3 2.1 Problemformulering . . . 3 2.2 Restriktioner . . . 3 2.3 Metoder . . . 4 2.3.1 Lokal modellering . . . 4 2.3.2 Linjär regression . . . 5 2.3.3 Lokal regression . . . 6 2.3.4 Direkt viktoptimering . . . 82.3.5 Multivariat direkt viktoptimering . . . 11
2.3.6 Regression med gaussiska processer . . . 13
2.4 Direkt viktoptimering med hjälp av Karush-Kuhn-Tucker-villkor . . . 14
3 Prestandamått och lokal skattning av inparameter 17 3.1 Globala prestandamått . . . 17
3.1.1 Absolutfel och genomsnittligt absolutfel . . . 17
3.1.2 Medelkvadratfelet och kvadratroten av medelkvadratfelet . 18 3.1.3 Korsvalidering . . . 18
3.1.4 Generaliserad korsvalidering . . . 19
3.2 Lokala prestandamått . . . 19
3.2.1 Lokal korsvalidering . . . 19
3.2.2 Lokal generaliserad korsvalidering . . . 20
3.2.3 Localized Akaike’s Information Criterion . . . 20
3.2.4 Localized Final Prediction Error . . . 21
3.3 Lokala prestandamått på direkt viktoptimering . . . 21
3.4 Lokal parameterskattning till direkt viktoptimering med hjälp av lokala prestandamått . . . 22
4 Implementering 23
4.1 Direkt viktoptimering . . . 23
4.1.1 Målfunktion . . . 24
4.1.2 Bivillkor . . . 24
4.2 Lokal parameterskattning till direkt viktoptimering . . . 25
5 Simuleringar 27 5.1 Utan brus . . . 27
5.1.1 Simulering 1: Utan brus . . . 27
5.2 Med brus . . . 29
5.2.1 Simulering 2: Med brus och varierad inparameter . . . 29
5.2.2 Simulering 3: Låg andraderivata . . . 30
5.2.3 Simulering 4: Hög andraderivata . . . 31
5.3 Monte Carlo-simuleringar . . . 31
5.3.1 Simulering 5, 6 och 7: Låg andraderivata och varierande brusnivå . . . 34
5.3.2 Simulering 8, 9 och 10: Hög andraderivata och varierande brusnivå . . . 36
5.3.3 Simulering 11, 12 och 13: Mycket hög andraderivata och va-rierande brusnivå . . . 36
5.4 Globalt skattad inparameter . . . 38
5.5 Val av den maximala inparametern λmax . . . 40
5.6 Finjusterad korsvalidering . . . 42
5.7 Prestanda . . . 43
5.8 Medelvärdesbildning . . . 44
5.9 Regression med gaussiska processer . . . 46
5.10 Ekvidistanta mätpunkter . . . 46
5.11 Multivariata funktioner . . . 46
6 Slutsatser 51 A Kod 55 A.1 Filen dwomultivariat.m . . . 55
Kapitel 1
Inledning
1.1
Bakgrund
Ämnet systemidentifiering är ett stort forskningsområde som under den senaste tiden växt snabbt. Systemidentifiering kan sägas handla om att ta fram matema-tiska modeller av system utifrån mätta data. Ett system kan här vara i stort sett vad som helst som har ett utfall, eller som vi kallar det en utsignal, som vi kan mäta. Denna utsignal beror på en eller flera faktorer där vi kan kontrollera ibland inga och ibland flera av dem. De faktorer som vi kan kontrollera kallar vi insignaler till systemet och de faktorer som vi inte kan kontrollera kallar vi störningar. Stör-ningarna delas ofta upp i mätbara och omätbara. Det är konsten att ta reda på hur utsignalen från systemet påverkas av insignalerna och/eller störningarna till systemet som kallas systemidentifiering, se t.ex. Ljung (1987) för mer information om systemidentifiering. De system där störningarna utgörs av additivt vitt brus är den typen av system som kommer att behandlas här. Kortfattat kan man säga att det som här är av intresse är att beskriva relationen mellan mätta insignaler och utsignaler så bra som möjligt.
Det finns många metoder för systemidentifiering men den metod på vilken tyngdpunkten i denna uppsats kommer att ligga är metoden från Roll (2003) som kallas direkt viktoptimering (Direct Weight Optimization, DWO). Problemet med denna metod är att den har en inparameter som innebär att man måste veta något om systemet i förväg vilket inte är önskvärt. Det finns dock metoder för att skatta denna inparameter och med idén från avsnitt 7.2 i Roll (2003) kan man få en lokal skattning av inparametern. Denna idé att lokalt skatta inparametern till DWO sammanfaller delvis med Model-on-Demand-tankesättet från Stenman (1999) där ett antal sätt att göra en lokal skattning föreslås. Metoder där man skattar en inparameter lokalt brukar man kalla adaptiva (ofta kallas skattningen Adaptive Bandwidth Selection). Dessa metoder har till fördel att de kan anpassa sig själva till till exempel funktioner med varierande första- eller andraderivata.
Att basera en adaptiv metod på DWO som i avsnitt 7.2 i Roll (2003) är till stora delar outforskat vilket är grunden till denna uppsats.
1.2
Syfte
Idén från Roll (2003) om lokalt skattad inparameter till DWO ger ett antal metoder att välja mellan och syftet med denna uppsats är således att ta reda på vilken av dessa metoder som ger det bästa resultatet. Syftet är också att undersöka om någon av dessa metoder alltid är bäst eller om det varierar från fall till fall och om det är någon idé att använda någon av dessa metoder eller om det är bättre att använda någon annan metod som inte bygger på DWO.
1.3
Disposition
Kapitel 2 och 3 innehåller bakgrundsteori för delarna av ämnet systemidentifiering som uppsatsen ska behandla. Kapitel 2 handlar främst om olika metoder som används vid systemidentifiering. Kapitel 3 behandlar olika sätt att mäta hur bra dessa metoder är och hur man kan använda detta till att skatta inparametern till DWO lokalt. I kapitel 4 finns information om hur man kan implementera DWO. Kapitel 5 är det huvudsakliga resultatkapitlet och innehåller jämförelser i form av simuleringar mellan de olika metoderna för att skatta inparamtern till DWO och några andra systemidentifieringsmetoder. Kapitel 6 innehåller slutsatser.
Kapitel 2
Systemidentifiering
I detta kapitel kommer vi fördjupa oss i ämnet systemidentifiering. Eftersom det vi i slutändan vill undersöka är om man kan förbättra eller förenkla resultat från Roll (2003) är det naturligt att stora delar av detta kapitel kommer vara likt delar ur just den avhandlingen.
2.1
Problemformulering
Här kommer vi behandla system där vi har mätt signalerna ϕ(k) och y(k). Vi kommer att kalla ϕ(k) för insignaler (kallas även regressionsvektor) och y(k) för utsignaler. Notera att ϕ(k) är insignaler till vår modell och de behöver alltså inte vara insignaler till systemet. Vi sätter ihop ϕ(k) och y(k) till insignal-utsignalpar {(ϕ(k), y(k))}N
k=1där N är antalet mätta värden.
Vi vill ställa upp en relation mellan dessa insignal-utsignalpar och vi vet att det kommer finnas brus från t.ex. mätutrustning som bör vara med i modellen. Här antar vi att allt brus är additivt. Ett sätt att ställa upp denna relation är genom modellen
y(k) = f (ϕ(k)) + e(k), (2.1)
där e(k) är individuellt oberoende likafördelade stokastiska variabler oberoende av regressionsvekorn ϕ(k) med E[e(k)] = 0 och E[e2(k)] = σ2. Här står E[ . ] för det statistiska väntevärdet. Funktionen f : Rn → R antas okänd och det är den som vi är intresserade av att skatta så bra som möjligt.
2.2
Restriktioner
Låt oss här betrakta det enklare fallet av (2.1) när f : R → R. Vi vill skatta f så bra som möjligt i punkten ϕ0 och betecknar denna skattning ˆf (ϕ0). För att våra mätvärden ska säga någonting om f måste vi sätta krav på hur snabbt f får ändra sig. Utan restriktion kan vi inte säga mycket om f , t.ex. om vi har mätt insignal-utsignalparen (−0.01, 0) och (0.01, 0) kan mycket väl f (0) vara 100.
Ett sätt att begränsa f är att anta att den är kontinuerligt deriverbar och att derivatan uppfyller ett så kallat lipschitzvillkor (se t.ex. Kreyszig (1989)), alltså att det finns en konstant L sådan att
|f0(ϕ(1)) − f0(ϕ(2))| ≤ L|ϕ(1) − ϕ(2)|, ∀ϕ(1), ϕ(2) ∈ R. (2.2) Vi kallar den klass av funktioner som uppfyller (2.2) F2(L). Att f ∈ F2(L) innebär således att f :s derivata bara får ändra sig ändligt fort. Detta medför att om vi har mätt (−0.01, 0) och (0.01, 0) och L är relativt liten vet vi att f (0) ligger nära 0.
2.3
Metoder
Det finns många sätt att hitta en skattning av f (ϕ0). I och med att vi ställde ett lipschitzkrav på f0 känns det naturligt att basera vår skattning på mätningar y(k) som ligger nära ϕ0.
2.3.1
Lokal modellering
Lokal modellering kallas det när man skattar f (ϕ0) genom att använda mätvär-den, y(k), som motsvarar värden som ligger nära ϕ0. En enkel metod, som kallas närmastegrannemetoden (nearest neighbor, se t.ex. Stenman (1999)), är att ta det värde på y(k) som motsvarar det ϕ(k) som ligger närmast ϕ0 som vårt ˆf (ϕ0). Uttryckt i formler använder man alltså
ˆ f (ϕ0) = y(m), (2.3) där m kommer från m = arg min m=1,...,N |ϕ(m) − ϕ0|,
där N som innan är totala antalet mätvärden. Denna metod tar bara med ett mätvärde och kommer därför vara mycket känslig för brus.
En utbyggnad av närmastegrannemetoden är att ta med fler värden och ta deras medelvärde. Denna metod kallas då k-närmastegrannemetoden (k-nearest neighbor, kNN, se t.ex. Stenman (1999)). För två stycken värden (k = 2) blir det
ˆ f (ϕ0) = y(m1) + y(m2) 2 , (2.4) där m1och m2kommer från m1= arg min m=1,...,N |ϕ(m) − ϕ0| (2.5) m2= arg min m=1,...,m1−1,m1+1,...,N |ϕ(m) − ϕ0|. (2.6)
Generellt, om kNN-metoden används blir skattningen
ˆ f (ϕ0) = k X i=1 1 ky(mi), (2.7)
där vi får mi genom att använda (2.5) och (2.6) och sedan fortsätta på samma
sätt.
Bruskänsligheten hos närmastegrannemetoden är dämpad hos kNN-metoden men kNN-metoden har istället problemet att den har svårt att hinna med snabba ändringar i f . Vi ser också att om k = N är inte detta en lokal metod längre utan en global.
Det är lätt att se att (2.7) kan skrivas som en viktad summa
ˆ f (ϕ0) = N X k=1 wky(k), (2.8)
med vikterna wk = 1/k då k = m1, . . . , mM och wk = 0 annars. Även (2.3) och
(2.4) kan skrivas på detta sätt med k = 1 respektive k = 2. Vi har alltså hittat ett sätt att skriva dessa metoder som en viktad summa av y-komponenterna.
Ett nästa steg är att använda en lite mer avancerad metod från statistiken.
2.3.2
Linjär regression
Tanken med linjär regression (se t.ex. Blom (1989)) är att ansätta en funktion som är linjär i ett antal okända parametrar mellan in- och utsignalen och sedan anpassa denna så bra som möjligt till mätta värden. Själva anpassningen utgörs oftast av en minstakvadratskattning.
En vanlig funktion att ansätta är ett polynom, t.ex. en rak linje eller en parabel. Vi kommer här anpassa ett polynom av grad p till våra mätta värden och sedan ur det bestämma ˆf (ϕ0). Vi vill alltså anpassa parametrarna i
P (ϕ) = β0+ β1ϕ + β2ϕ2+ · · · + βpϕp, (2.9)
så bra som möjligt till vårt insignal-utsignalpar {(ϕ(k), y(k))}N
k=1, där y(k) fås
som
y(k) = P (ϕ(k)) + e(k), (2.10)
där e(k) är felet mellan mätningen och skattningen. Vi inför Y = (y(1), . . . , y(N ))T, β = (β0, . . . , β
p)T, E = (e(1), . . . , e(N ))T och
Φ = 1 ϕ(1) · · · ϕ(1)p .. . ... ... 1 ϕ(N ) · · · ϕ(N )p .
Med detta, (2.9) och (2.10) ser vi att
Y = y(1) .. . y(N ) = β0+ β1ϕ(1) + · · · + βpϕ(1)p .. . β0+ β1ϕ(N ) + · · · + βpϕ(N )p + e(1) .. . e(N ) = 1 ϕ(1) · · · ϕ(1)p .. . ... ... 1 ϕ(N ) · · · ϕ(N )p β0 .. . βp + e(1) .. . e(N ) = Φβ + E, (2.11)
där N antas vara större än p. Vi är intresserade av att bestämma det β som minimerar felet E och det är därför naturligt att använda minstakvadratmetoden. Minstakvadratlösningen betecknar vi med ˆβ och vi får, så länge (ΦTΦ)−1existerar,
ˆ β = (ΦTΦ)−1ΦTY. (2.12) Vidare sätter vi Ψ(ϕ) = 1 ϕ .. . ϕp .
Polynomet (2.9) med anpassade parametrar fås då som
ˆ
P (ϕ) = Ψ(ϕ)Tβ.ˆ
Detta polynom i punkten ϕ0 är vårt sökta ˆf (ϕ0). Vi får alltså slutligen
ˆ f (ϕ0) = ΨT(ϕ0)(ΦTΦ)−1 N X k=1 Ψ(ϕ(k))y(k) = N X k=1 wky(k), där wk= ΨT(ϕ0)(ΦTΦ)−1Ψ(ϕ(k)).
Vi ser därmed att även denna lite mer avancerade metod slutar i en viktad summa av y-komponenterna likt (2.8). Den enda skillnaden är att vikterna är olika.
2.3.3
Lokal regression
I avsnitt 2.3.2 använde vi minstakvadratmetoden för att anpassa ett polynom till våra mätvärden. Problemet med det angreppssättet är att man värderar värden långt från ϕ0 lika mycket som värden nära ϕ0. Detta försöker man mildra i den lokala varianten av linjär regression genom att införa en viktfunktion, K(ϕ), som man då kan välja så att värden nära ϕ0 ges större prioritet.
Vi vet från minstakvadratmetoden att (2.12) kan skrivas
ˆ β = arg min β N X k=1 y(k) − p X i=0 βiϕi(k) !2 .
Om vi viktar denna med vårt K(ϕ) och skattar runt ϕ0kallas det lokal regression (Local Polynomial Regression, LPR, eller Local Polynomial Modelling, se t.ex. Fan och Gijbels (1996)). Parameterskattningen betecknas ˆβ och parametrarna
anges av β = (β0, . . . , βp)T. Eftersom vi vill skatta runt ϕ0 blir ϕ(k) ersatt med ˜
ϕ(k) = ϕ(k) − ϕ0. Alltså kan skattningen skrivas
ˆ β = arg min β N X k=1 K( ˜ϕ(k)) y(k) − p X i=0 βiϕ˜i(k) !2 . (2.13)
Om vi inför Y = (y(1), . . . , y(N ))T, Φ = 1 ϕ(1)˜ · · · ϕ˜p(1) .. . ... ... 1 ϕ(N )˜ · · · ϕ˜p(N ) (2.14) och ¯ K = K( ˜ϕ(1)) 0 · · · 0 0 K( ˜ϕ(2)) . . . 0 .. . ... . .. ... 0 0 . . . K( ˜ϕ(N ))
kan vi skriva (2.13) på matrisform genom
ˆ
β = arg min
β
(Y − Φβ)TK(Y − Φβ),¯
som är en lokal minstakvadratlösning (se t.ex. Hackman (2001)) som, så länge (ΦTKΦ)¯ −1 existerar, kan skrivas
ˆ
β = (ΦTKΦ)¯ −1ΦTKY.¯
Vi är intresserade av ˆf (ϕ0) och vet att med ˜ϕ = ϕ − ϕ0 är, likt (2.9),
ˆ
P (ϕ) = ˆβ0+ ˆβ1ϕ + ˆ˜ β2ϕ˜2+ · · · + ˆβpϕ˜p,
som i punkten ϕ0 är just detta ˆf (ϕ0). Eftersom ˜ϕ = 0 när ϕ = ϕ0 får vi
ˆ
f (ϕ0) = ˆβ0.
Vi kan alltså få vårt sökta ˆf (ϕ0) genom
ˆ
f (ϕ0) = ˆβ0= eT1β = eˆ
T
1(Φ
TKΦ)¯ −1ΦTKY,¯ (2.15)
där e1är enhetsvektorn (1, 0, . . . , 0)T. Vidare vill vi skriva (2.15) på formen
ˆ f (ϕ0) = N X k=1 wky(k) = wTY . Detta ger ˆ f (ϕ0) = wTY = eT1(Φ TKΦ)¯ −1ΦTKY =¯ eT 1(Φ TKΦ)¯ −1ΦTK¯TT Y = K¯TΦ((ΦTKΦ)¯ T)−1e1 T Y = KΦ(Φ¯ TKΦ)¯ −1e1 T Y, (2.16)
där vi i den sista likheten utnyttjat att både ¯K och ΦTKΦ är symmetriska. Vi¯ har alltså
som ger
wk= eTkKΦ(Φ¯
TKΦ)¯ −1e
1.
Vi kan alltså slutligen skriva vår skattade funktion som
ˆ f (ϕ0) = N X k=1 wky(k), där wk= eTkKΦ(Φ¯ TKΦ)¯ −1e 1.
Alltså kan resultatet även med denna metod skrivas som en viktad summa av y-komponenterna. Att många metoder slutar i en viktad summa kan ses som en motivering till nästa metod.
2.3.4
Direkt viktoptimering
I direkt viktoptimering (Direct Weight Optimization, DWO, Roll (2003)), skattar man ˆf (ϕ0) genom att starta med en viktad summa
ˆ f (ϕ0) = N X k=1 wky(k) (2.18)
och sedan direkt, utan mellansteg, optimera fram vikterna, wk, genom att
mini-mera WMSE (Worst-case Mean Square Error) med avseende på wk.
För att få WMSE definierar vi först MSE (Mean Square Error) som
M SE ˆf (ϕ0) , E[ ˆf (ϕ0) − f (ϕ0) 2
|{ϕk}Nk=1]. (2.19)
Vi får sedan WMSE genom att ta supremum av MSE över en funktionsklass F , det vill säga vi definierar WMSE som
W M SE ˆf (ϕ0), F , sup
f ∈F
E[ ˆf (ϕ0) − f (ϕ0) 2
|{ϕk}Nk=1]. (2.20)
(2.18) och y(k) enligt (2.1). På samma sätt som i Roll (2003) får vi W M SE ˆf (ϕ0), F2(L) = sup f ∈F2(L) E[ N X k=1 wky(k) − f (ϕ0) !2 |{ϕk}Nk=1] = sup f ∈F2(L) E[ N X k=1 wk f (ϕ(k)) + e(k)− f (ϕ0) !2 |{ϕk}Nk=1] = sup f ∈F2(L) N X k=1 wkf (ϕ(k)) − f (ϕ0) !2 + σ2 N X k=1 wk2 (2.21) = sup f ∈F2(L) f (ϕ0) N X k=1 wk− 1 ! + f0(ϕ0) N X k=1 wkϕ(k)˜ + N X k=1 wk f (ϕ(k)) − f (ϕ0) − f0(ϕ0) ˜ϕ(k) !2 + σ2 N X k=1 w2k,
där vi i den sista likheten har lagt till och dragit bort
f (ϕ0) N X k=1 wk+ f0(ϕ0) N X k=1 wkϕ(k).˜
Detta har gjorts för att vi ska kunna använda följande sats (lemma 4.1.12 i Dennis och Schnabel (1996)), där ∇f står för gradienten till f och k · k2 står för den Euklidiska normen. Satsen anges här i vektorform för att kunna användas i avsnitt 2.3.5.
Sats 2.1 Låt f : Rn→ R vara en kontinuerligt deriverbar funktion vars derivata
uppfyller Lipschitzvillkoret k∇f (ϕ(k)) − ∇f (ϕ0)k2≤ Lk ˜ϕ(k)k2, ∀ϕ(k), ϕ0∈ Rn. (2.22) Då är |f (ϕ(k)) − f (ϕ0) − ∇Tf (ϕ0) ˜ϕ(k)| ≤L 2k ˜ϕ(k)k 2 2. (2.23)
Bevis Vi inför först g(t) = ϕ0+ t ˜ϕ och A = ∇f . Vi noterar sedan att
ϕ0+ ˜ϕ Z ϕ0 A · dr = ϕ0+ ˜ϕ Z ϕ0 ∇T f (r)dr = 1 Z 0 ∇T f (g(t))dg(t) dt dt = 1 Z 0 ∇Tf (g(t)) ˜ϕdt = 1 Z 0 ∇Tf (ϕ0+ t ˜ϕ) ˜ϕdt, (2.24)
där ”·” står för skalärprodukten. Med hjälp av (2.24) och vektoranalysens teori om lin-jeintegral och gradient (se t.ex. Ramgard (2000)) får vi
1 Z 0 ∇T f (ϕ0+ t ˜ϕ) ˜ϕdt = ϕ0+ ˜ϕ Z ϕ0 A · dr = f (ϕ0+ ˜ϕ(k)) − f (ϕ0). Med hjälp av detta får vi |f (ϕ(k)) − f (ϕ0) − ∇Tf (ϕ0) ˜ϕ(k)| = |f (ϕ0+ ˜ϕ(k)) − f (ϕ0) − ∇Tf (ϕ0) ˜ϕ(k)| = 1 Z 0 ∇Tf (ϕ0+ t ˜ϕ(k)) − ∇Tf (ϕ0)ϕ(k) dt˜ ≤ 1 Z 0 k∇f (ϕ0+ t ˜ϕ(k)) − ∇f (ϕ0)k2k ˜ϕ(k)k2dt ≤ 1 Z 0 Ltk ˜ϕ(k)k22dt = L 2k ˜ϕ(k)k 2 2. (2.25)
I Roll (2003) finns ett lemma, Lemma A.3, som är en mer allmän version av sats 2.1 med ett annat bevis.
I vårt fall, när f : R → R, blir (2.22) |f0(ϕ(k)) − f0(ϕ0)| ≤ L| ˜ϕ(k)|, ∀ϕ(k), ϕ0∈ R och (2.23) blir |f (ϕ(k)) − f (ϕ0) − f0(ϕ0) ˜ϕ(k)| ≤ L 2ϕ(k)˜ 2.
Enligt avsnitt 2.2 uppfyller f kraven och vi kan således använda satsen.
Innan vi använder satsen ser vi i (2.21) att eftersom vi tar supremum över funktionsklassen F2(L) kan f (ϕ0) och f0(ϕ0) vara godtyckligt stora. Vi måste därför se till att N X k=1 wk = 1 (2.26a) N X k=1 wkϕ(k) = 0.˜ (2.26b)
om ˆf (ϕ0) i (2.18). Detta, (2.21) och den univariata varianten av sats 2.1 ger W M SE ˆf (ϕ0), F2(L) = sup f ∈F2(L) N X k=1 wk f (ϕ(k)) − f (ϕ0) − f0(ϕ0) ˜ϕ(k) !2 + σ2 N X k=1 w2k ≤ sup f ∈F2(L) N X k=1 |wk| f (ϕ(k)) − f (ϕ0) − f 0(ϕ0) ˜ϕ(k) !2 + σ2 N X k=1 w2k ≤L 2 4 N X k=1 ˜ ϕ(k)2|wk| !2 + σ2 N X k=1 wk2. (2.27)
Vi har nu kommit fram till någonting vi kan minimera och minimerar alltså en övre gräns till WMSE. Tillsammans med (2.26) kan vi skriva minimeringen av (2.27) som ett konvext minimeringsproblem med bivillkor,
min w L2 4 N X k=1 ˜ ϕ2(k)|wk| !2 + σ2 N X k=1 wk2 då N X k=1 wk= 1 N X k=1 wkϕ(k) = 0,˜ (2.28)
vars lösning vi skriver w∗k. Vi får sedan vårt sökta ˆf (ϕ0) genom att sätta in w∗k i (2.18).
2.3.5
Multivariat direkt viktoptimering
Att anta att funktionen f är univariat räcker ofta inte eftersom system ofta har mer än en insignal. I detta fall vill man skapa en DWO-metod som försöker skatta en multivariat funktion, alltså en funktion f : Rn→ R. Detta kan i sin tur utvidgas till system med flera utsignaler genom att man använder DWO-metoden en gång för varje utsignal.
Vi kommer här att anta att f har samma restriktioner som anges i avsnitt 2.2 men med skillnaden att f : Rn → R. Detta medför att istället för kravet (2.2) har
vi
k∇f (ϕ(k)) − ∇f (ϕ0)k2≤ Lk ˜ϕ(k)k2, ∀ϕ(k), ϕ0∈ Rn.
Man ser lätt att många saker från avsnitt 2.3.4 kommer vara mycket lika i det multivariata fallet. Skillnaderna i kalkylen (2.21) blir att f0 ersätts med ∇Tf och
att alla ϕ är vektorer. Alltså W M SE ˆf (ϕ0), F2(L) = sup f ∈F2(L) f (ϕ0) N X k=1 wk− 1 ! + ∇Tf (ϕ0) N X k=1 wkϕ(k)˜ + N X k=1 wk f (ϕ(k)) − f (ϕ0) − ∇Tf (ϕ0) ˜ϕ(k) !2 + σ2 N X k=1 wk2. (2.29)
Utifrån detta kan vi se att kraven (2.26) blir
N X k=1 wk= 1 (2.30a) N X k=1 wkϕ˜i(k) = 0, i = 1, . . . , n. (2.30b)
Dessa krav, sats 2.1 och (2.29) ger likt (2.27)
W M SE ˆf (ϕ0), F2(L) = sup f ∈F2(L) N X k=1 wk f (ϕ(k)) − f (ϕ0) − ∇Tf (ϕ0) ˜ϕ(k) !2 + σ2 N X k=1 wk2 ≤L 2 4 N X k=1 k ˜ϕ(k)k22|wk| !2 + σ2 N X k=1 wk2.
Detta ger oss minimeringsproblemet
min w L2 4 N X k=1 k ˜ϕ(k)k22|wk| !2 + σ2 N X k=1 wk2 då N X k=1 wk= 1 N X k=1 wkϕ˜i(k) = 0, i = 1, . . . , n, (2.31)
skriva om detta till min w,s L2 4 N X k=1 k ˜ϕ(k)k22sk !2 + σ2 N X k=1 s2k (2.32a) då sk≥ wk sk≥ −wk (2.32b) N X k=1 wk = 1 N X k=1 wkϕ˜i(k) = 0, i = 1, . . . , n. (2.32c)
Att detta är samma som (2.31) kan man se genom att först notera att de två första bivillkoren i (2.32) tillsammans leder till sk≥ |wk|. Att både k ˜ϕ(k)k22och σ2alltid är positiva innebär att om sk är större än wk kan vi alltid minska sk utan att
bivillkoren bryts vilket leder till att målfunktionen minskar. När vi har minskat
sk så mycket som bivillkoren tillåter är sk= |wk| och vi ser därmed att (2.31) och
(2.32) är ekvivalenta.
Vi kan se att om vi delar målfunktionen i (2.32) med σ2beror problemet bara av en konstant, L/σ. Denna kvot är alltså det enda vi behöver veta om systemet och om noggrannheten i mätningarna när vi försöker skatta en funktion med hjälp av DWO. Vi inför därför λ = L/σ.
2.3.6
Regression med gaussiska processer
Regression med gaussiska processer (GP, se t.ex. Rasmussen och Williams (2006)), kallas även kriging (Krige (1951)) är en modern regressionsmetod som man får genom att man antar att mätvärdena kommer från en så kallad gaussisk process (se t.ex. Yates och Goodman (1999)). En egenskap som gaussiska processer har är att de kan beskrivas fullständigt av sin väntevärdesfunktion
µ(ϕ) = E[f (ϕ)]
och kovariansfunktion
Cov[f (ϕ(i)), f (ϕ(j))] = E[(f (ϕ(i)) − µ(ϕ(i)))(f (ϕ(j)) − µ(ϕ(j)))]
= C(ϕ(i), ϕ(j)) (2.33)
Här antas alltså att alla f (ϕ(k)), k = 1, . . . , N kommer från en gaussisk process. Man brukar också ofta anta att µ(ϕ) är noll.
Vi väljer att använda samma modell som i avsnitt 2.1 med det extra kravet att bruset måste vara normalfördelat. Vi har alltså brusvarians σ2. Om vi då inför
CN = C(ϕ(1), ϕ(1)) C(ϕ(1), ϕ(2)) . . . C(ϕ(1), ϕ(N )) C(ϕ(2), ϕ(1)) C(ϕ(2), ϕ(2)) . . . C(ϕ(2), ϕ(N )) .. . ... . .. ... C(ϕ(N ), ϕ(1)) C(ϕ(N ), ϕ(2)) . . . C(ϕ(N ), ϕ(N )) ,
visar det sig (se t.ex. Rasmussen och Williams (2006)) att man kan härleda funk-tionsskattningen ˆ f (ϕ0) = N X k=1 vkC(ϕ(k), ϕ0), (2.34) där v = v1 .. . vN = (CN + σ 2I)−1Y. Vidare inför vi Q = (CN + σ2I)−1
och kan då skriva (2.34) som
ˆ f (ϕ0) = N X k=1 wkyk, där wk= N X i=1 Q(i, k)C(ϕ(i), ϕ0).
GP kan alltså även den ses som en viktad summa av y-komponenterna. Vikterna är beroende av kovariansfunktionen som inte är given från början. Man måste alltså välja en kovariansfunktion i förväg och detta val kan man läsa mer om i t.ex. Rasmussen och Williams (2006). Att GP är en modern regressionsmetod som bygger på en viktad summa gör att den blir intressant att jämföra DWO med.
2.4
Direkt viktoptimering med hjälp av
Karush-Kuhn-Tucker-villkor
Karush-Kuhn-Tucker-villkor (KKT-villkor) är optimalitetsvillkor som framför allt används till problem med ickelinjär målfunktion och ickelinjära bivillkor (se t.ex. Lundgren m.fl. (2003)). I avsnitt 3.2 i Roll (2003) visas att om vi använder KKT på DWO och numrerar om våra vikter w∗t från DWO så att alla w∗k 6= 0 då
k = 1, . . . , M och inför rk= 1, w∗k> 0 −1, w∗k< 0 , k = 1, . . . , M, ¯ ˜ ϕ = ˜ ϕT(1) .. . ˜ ϕT(M ) , z = r1k ˜ϕ(1)k22 .. . rMk ˜ϕ(M )k22 , w ∗= w∗1 .. . w∗ M och Φ = (1M ϕ),¯˜ G = λ2 4 zz T + I, (2.35)
Karush-Kuhn-Tucker-villkor 15
kan vi skriva
w∗= G−1Φ(ΦTG−1Φ)−1e1.
Om vi jämför detta med (2.17) och vårt Φ med (2.14) ser vi att vikterna från DWO som inte är noll kan ses som vikterna från lokal regression av grad ett, där
¯
K i (2.17) blivit ersatt av G−1 från (2.35). Tyvärr kan vi inte använda detta till att räkna ut w∗ direkt eftersom vi inte kan veta rk eller vilka vikter som blir noll
Kapitel 3
Prestandamått och lokal
skattning av inparameter till
direkt viktoptimering
I kapitel 2 presenterades några olika metoder för att skatta en funktion och den vi är mest intresserade av här är DWO-metoden. När man väl har gjort en skattning vill man också ha något sätt att veta hur bra skattningen är. I detta kapitel kommer därför ett antal prestandamått att presenteras. Vi kommer också att undersöka hur man kan använda dessa mått för att lokalt skatta ett λ som kan användas i DWO-metoden.
I alla prestandamått i detta kapitel anges hur bra skattningen har varit med en siffra där lägre betyder bättre. De fyra första måtten är enkla standardmått som används i många sammanhang.
3.1
Globala prestandamått
Med globala prestandamått menas här mått som mäter hur bra hela vår skattade funktion ˆf är. Vi börjar med de två enklaste måtten som är mycket lika varandra.
3.1.1
Absolutfel och genomsnittligt absolutfel
Absolutfel (Absolute Error, AE) är helt enkelt summan av felen i skattningarna man får om man skattar f i punkterna ϕ(k), k = 1, . . . , N . Vi får
AE( ˆf ) = N X k=1 f (ϕ(k)) − ˆf (ϕ(k)) . (3.1)
Det är ofta mer intressant att veta hur stort felet är i medel med avseende på
N . Detta får man i genomsnittligt absolutfel (Mean Absolute Error, MAE), som
då blir M AE( ˆf ) = 1 N N X k=1 f (ϕ(k)) − ˆf (ϕ(k)) . (3.2)
3.1.2
Medelkvadratfelet och kvadratroten av
medelkvadrat-felet
Medelkvadratfelet (Mean Square Error, MSE), som vi först såg i avsnitt 2.3.4, är likt MAE men där felet är kvadrerat för att se till att stora fel får större vikt. MSE får vi här genom M SE( ˆf ) = 1 N N X k=1 f (ϕ(k)) − ˆf (ϕ(k)) 2 . (3.3)
Kvadratroten av medelkvadratfelet (Root Mean Square Error, RMSE) är helt enkelt kvadratroten av MSE. Detta för att det ska få samma enhet som f , vilket ibland kan vara av intresse. Alltså gäller det att
RM SE( ˆf ) = v u u t PN k=1 f (ϕ(k)) − ˆf (ϕ(k)) 2 N . (3.4)
3.1.3
Korsvalidering
I (3.1), (3.2), (3.3) och (3.4) används det sanna f som ju är det vi försöker skatta och således inte alltid kan veta. Man kan därför försöka skatta f (ϕ(k)) med det mätta y(k), vilket till exempel i MSE ger
M SEy( ˆf ) = 1 N N X k=1 y(k) − ˆf (ϕ(k)) 2 , (3.5)
där M SEy står för det skattade MSE som vi får med hjälp av y(k). Ett problem
kan dock vara att ˆf (ϕ(k)) är en viktad summa av y(k) i DWO-fallet. Om vi då
väljer vikter så att ˆf (ϕ(k)) = y(k) blir M SEy= 0. Det betyder att vi alltid kan få
en perfekt skattning i M SEy-mening, men att skattningen då blir mycket känslig
för brus. Detta är alltså inte ett bra sätt att skatta MSE.
Man kan gå runt detta problem genom att använda så kallad korsvalidering (Leave-One-Out Cross-Validation, CV, se t.ex. Stenman (1999) sid. 38) som går ut på att man ersätter ˆf (ϕ(k)) med den skattning man skulle få om man utelämnar
mätning nummer k. Vi kallar denna skattning ˆf−k(ϕ(k)) och får
CV ( ˆf ) = 1 N N X k=1 y(k) − ˆf−k(ϕ(k)) 2 . (3.6)
Om man inte vill räkna ut ˆf−k(ϕ(k)) explicit kan man använda (se t.ex. Stenman (1999) sid. 39) CV ( ˆf ) = 1 N N X k=1 y(k) − ˆf (ϕ(k)) 1 − infl(ϕ(k)) !2 , (3.7)
där infl(ϕ(k)) är den så kallade Influence Function (se t.ex. Stenman (1999) sid. 35) som om man har lokal regression och beteckningarna från avsnitt 2.3.3 är
infl(ϕ0) = eT1(Φ
TKΦ)¯ −1e
1K(0). (3.8)
3.1.4
Generaliserad korsvalidering
I den så kallade generaliserade korsvalideringen (Generalized Cross-Validation, GCV, se t.ex. Stenman (1999) sid. 39), byter man infl(ϕ(k)) i (3.7) mot tr(H)/N , där tr(·) står för spåret och H = wT(ϕ(1)) .. . wT(ϕ(N ) , (3.9)
där w(ϕ) står för vikterna man får om man skattar funktionen f i punkten ϕ. Man kan se tr(H)/N som medelvärdet av alla infl(ϕ(k)) (se t.ex. Stenman (1999) sid. 39). GCV blir alltså
GCV ( ˆf ) = N N X k=1 y(k) − ˆf (ϕ(k)) N − tr(H) !2 . (3.10)
3.2
Lokala prestandamått
Med lokala prestandamått menas här mått som mäter hur bra vår skattning är i en viss punkt, alltså hur bra ˆf (ϕ0) är. Vi kommer här se hur CV och GCV ser ut i sina lokala tappningar och införa två nya mått.
3.2.1
Lokal korsvalidering
Den lokala versionen av CV som vi kallar lokal korsvalidering (Localized Cross-Validation, LCV, se t.ex. Stenman (1999) sid. 57) förutsätter att ˆf kommer från
en lokal regression som i avsnitt 2.3.3. Med beteckningarna från avsnitt 2.3.3 får vi LCV ( ˆf (ϕ0)) = PN k=1K( ˜ϕ(k)) y(k) − ¯f ϕ0 −k(ϕ(k)) 2 PN k=1K( ˜ϕ(k)) , (3.11)
där ¯fϕ0
−k(ϕ) står för det skattade f vi får i punkten ϕ om vi utelämnar mätning nummer k i en lokal regression. Om vi inför
F1:N = ¯ fϕ0 −k(ϕ(1)) .. . ¯ fϕ0 −k(ϕ(N )) , (3.12)
kan vi skriva (3.11) som
LCV ( ˆf (ϕ0)) =
(Y − F1:N)TK(Y − F¯ 1:N)
tr( ¯K) . (3.13)
I avsnitt 2.4 såg vi att DWO är en lokal regression och vi kan alltså använda LCV för att få ett lokalt prestandamått för DWO.
3.2.2
Lokal generaliserad korsvalidering
Lokal generaliserad korsvalidering (Localized Generalized Cross-Validation, LGCV, se t.ex. Stenman (1999) sid. 58) är en lokal variant av GCV. Notera att uttrycken skiljer sig åt men detta är troligtvis på grund av ett tryckfel i Stenman (1999). LGCV kan skrivas som
LGCV ( ˆf (ϕ0)) = tr( ¯K) PN k=1K( ˜ϕ(k)) (y(k) − ¯mk) 2 (tr( ¯K) − tr (ΦTKΦ)¯ −1ΦTK¯2Φ)2, (3.14) där ¯ mk= ekΦ(ΦTKΦ)¯ −1ΦTKY.¯ (3.15)
Om vi skriver detta på matrisform får vi
LGCV ( ˆf (ϕ0)) = tr( ¯K) (Y − ¯m)TK(Y − ¯¯ m) (tr( ¯K) − tr (ΦTKΦ)¯ −1ΦTK¯2Φ)2, (3.16) där ¯ m = Φ(ΦTKΦ)¯ −1ΦTKY.¯ (3.17)
Även LGCV kan användas på en funktionsskattning från en lokal regression vilket leder till att vi kan använda detta på DWO.
3.2.3
Localized Akaike’s Information Criterion
En lokaliserad variant av ett annat vanligt förekommande prestandamått är Loca-lized Akaike’s Information Criterion (LAIC, se t.ex. Stenman (1999) sid. 59) som i matrisform blir LAIC( ˆf (ϕ0)) = (Y − ¯m) TK(Y − ¯¯ m) tr( ¯K) exp α tr (ΦTKΦ)¯ −1ΦTK¯2Φ tr( ¯K) ! , (3.18)
3.2.4
Localized Final Prediction Error
En variant av LAIC är Localized (Akaike’s) Final Prediction Error (LFPE, se t.ex. Stenman (1999) sid. 59) som i matrisform blir
LF P E( ˆf (ϕ0)) =(Y − ¯m) TK(Y − ¯¯ m) tr( ¯K) × 2 tr( ¯K) + α tr (Φ TKΦ)¯ −1ΦTK¯2Φ 2 tr( ¯K) − α tr (ΦTKΦ)¯ −1ΦTK¯2Φ ! , (3.19)
där α, som i (3.18), är ett valfritt straff på varianstermen. Till skillnad från LAIC så saknas det instruktioner om hur man ska välja α i Stenman (1999).
Om vi tittar på (3.19) ser vi att ett för högt värde på α kan resultera i ne-gativa värden så det vore rimligt att välja α så att måttet garanterat är posi-tivt. Ett negativt värde på ett prestandamått är både mycket ovanligt och oin-tuitivt. Till exempel borde en negativ nämnare i (3.19) göra att en skattning långt från det mätta värdet ger ett till beloppet stort negativt LFPE-värde. Det innebär att om vi minimerar måttet borde vi få en så dålig anpassning som möj-ligt. Vidare analys ger att α ≤ 2 garanterar ett ickenegativt värde på måttet då 0 ≤ tr (ΦTKΦ)¯ −1ΦTK¯2Φ /tr( ¯K) ≤ 1 (se t.ex. Stenman (1999)). Om man dess-utom tänker på att α är ett straff på varianstermen borde det enda rimliga vara att
α är ickenegativ. Det är alltså rimligt att α ska väljas som 0 ≤ α ≤ 2. Det antyds
dock på sid. 65 i Stenman (1999) att α ≥ 2 ska väljas. Detta, som vi kommer att se senare, har också gett intressanta resultat trots att måttet ibland är negativt. Vi kommer därför i fortsättningen att använda α ≥ 0 i LFPE.
Även LAIC och LFPE kan användas på en lokal regression, alltså kan vi an-vända dem också på DWO.
3.3
Lokala prestandamått på direkt
viktoptime-ring
Vi vet alltså att vi kan använda LCV, LGCV, LAIC och LFPE för att få ett lokalt prestandamått på DWO. Vi såg också i avsnitt 2.4 att DWO kan ses som om det vore en lokal regression av grad ett, där ¯K i (2.17) har blivit ersatt av G−1 från (2.35). Det enda vi behöver göra är alltså att ersätta ¯K i (3.13), (3.16), (3.18)
beteckningarna från avsnitt 2.4 får vi LCV ( ˆf (ϕ0)) =(Y − F1:N) TG−1(Y − F1:N) tr(G−1) (3.20) LGCV ( ˆf (ϕ0)) =tr(G−1) (Y − ¯m)TG−1(Y − ¯m) (tr(G−1) − tr ((ΦTG−1Φ)−1ΦT(G−1)2Φ))2 (3.21) LAIC( ˆf (ϕ0)) =(Y − ¯m) TG−1(Y − ¯m) tr(G−1) × exp α tr (Φ TG−1Φ)−1ΦT(G−1)2Φ tr(G−1) ! (3.22) LF P E( ˆf (ϕ0)) = (Y − ¯m)TG−1(Y − ¯m) tr(G−1) × 2 tr(G −1) + α tr (ΦTG−1Φ)−1ΦT(G−1)2Φ 2 tr(G−1) − α tr ((ΦTG−1Φ)−1ΦT(G−1)2Φ) ! . (3.23)
3.4
Lokal parameterskattning till direkt
viktopti-mering med hjälp av lokala prestandamått
I avsnitt 2.3.5 såg vi att kvoten λ = L/σ var det enda vi behövde veta på förhand för att använda DWO. Vi är därför intresserade av att på något sätt välja ett λ utan att på förhand veta just detta om systemet och mätningarna.
Eftersom (3.20), (3.21), (3.22) och (3.23) kan användas för att få prestandamått på DWO och vi vet att de innehåller λ kan vi minimera dessa med avseende på λ och det erhållna ˆλ är alltså vårt sökta λ som kan användas i DWO-metoden. Om
vi vill använda LCV får vi alltså
ˆ
λ = arg min
λ
(LCV ( ˆf (ϕ0))), (3.24)
där LCV fås från (3.20). Om vi vill använda LGCV, LAIC eller LFPE fås ˆλ på
samma sätt men med LCV utbytt mot motsvarande prestandamått.
Vi har alltså slutligen kommit fram till något som går att använda för att lokalt skatta inparametern till DWO i skattningspunkten ϕ0. Detta är dock ännu inte löst analytiskt utan vi får förlita oss på simuleringar. Det vi är intresserade av att veta är hur bra denna metod fungerar och vilken av LCV, LGCV, LAIC eller LFPE som fungerar bäst. Detta kommer vi titta närmare på i kapitel 5.
Kapitel 4
Implementering
I avsnitt 3.3 kom vi fram till det sista vi behövde i teoriväg för att börja simu-lera. Det återstår dock några problem för att få det praktiskt genomförbart och någorlunda smidigt. Detta ska vi se hur vi kan lösa i detta kapitel.
4.1
Direkt viktoptimering
Vi såg i avsnitt 2.3.5 att det centrala i DWO-metoden är att lösa minimeringspro-blemet min w,s λ2 4 N X k=1 k ˜ϕ(k)k22sk !2 + N X k=1 s2k (4.1a) då sk ≥ wk sk ≥ −wk (4.1b) N X k=1 wk= 1 N X k=1 wkϕ˜i(k) = 0 i = 1, . . . , n. (4.1c)
Vi ser att i och med vår omskrivning i avsnitt 2.3.5 från (2.31) till (2.32), som är samma som (4.1) får vi ett så kallat Quadratic Problem (QP). Ett QP är ett minimeringsproblem som kan skrivas på formen
min x 1 2x TQx + cTx (4.2a) då Ax ≤ b (4.2b) Ex = d, (4.2c)
där x, c, b samt d är vektorer och Q, A samt E är matriser där Q är symmetrisk.
Att skriva om problemet på QP-form har till fördel att det finns effektiva lösare att tillgå för den typen av optimeringproblem. Vi kommer här att använda lösaren qpas från projektet QPC (Wills (2007)).
Det vi är ute efter är alltså x, Q, c, A, b, E och d så att vi kan skriva (4.1) på formen (4.2).
4.1.1
Målfunktion
Vi börjar med att sätta
x = s w , φ = λ 2 k ˜ϕ(1)k2 2 .. . k ˜ϕ(N )k22 , (4.3)
där s = (s1. . . sN)T och w = (w1. . . wN)T. Därefter sätter vi
Q = 2 φφT + IN 0N ×N 0N ×N 0N ×N , c = 02N ×1, (4.4)
där 0N ×N är en N × N matris med bara nollor och c är således en nollvektor med
2N element. Detta insatt i (4.2a) ger
1 2x TQx + cTx = 1 22 s T wT φφT+ I 0 N ×N 0N ×N 0N ×N s w + 0 =λ 2 4 s TφφTs + sTs =λ 2 4 N X k=1 k ˜ϕ(k)k22sk !2 + N X k=1 s2k.
Därmed ser vi att om vi väljer x, φ, Q och c enligt (4.3) och (4.4) kan målfunktionen (4.1a) skrivas som (4.2a). Vi ser att vi som önskat kan använda λ som enda obekant men om σ = 0 måste vi istället använda en målfunktion som i (2.32). Vi får då
φ =L 2 k ˜ϕ(1)k2 2 .. . k ˜ϕ(N )k2 2 , Q = 2 φφT 0 N ×N 0N ×N 0N ×N . (4.5)
4.1.2
Bivillkor
De två första bivillkoren, (4.1b), får vi genom att sätta
A = −IN IN −IN −IN , b = 02N ×1, (4.6)
där IN är identitetsmatrisen av storlek N . Bivillkoret (4.2b) ger då
Ax = −IN IN −IN −IN s w = IN(−s + w) IN(−s − w) ≤ 02N ×1= b.
Med andra ord sk ≥ wk och sk ≥ −wk. Därmed ser vi att (4.1b) kan skrivas som
(4.2b) med x, A och b enligt (4.3) och (4.6).
För att få de andra två bivillkoren, (4.1c), sätter vi ˜ϕ = ( ˜ϕ(1), . . . , ˜ϕ(N ))T,
d = (1, 0n×1)T och E = 01×N 11×N 0n×N ϕ˜T . (4.7) Då blir bivillkoren (4.2c) Ex = 01×N 11×N 0n×N ϕ˜T s w = PN k=1wk PN k=1wkϕ(k)˜ ! = 1 0n×1 = d.
Var för sig blir det
N X k=1 wk= 1 N X k=1 wkϕ˜i(k) = 0 i = 1, . . . , n,
det vill säga samma som i (4.1c). Vi kan alltså skriva (4.1c) som (4.2c) om vi väljer
x, E och d enligt (4.3) och (4.7).
Vi har alltså skrivit om (4.1) på formen (4.2) som vi var ute efter. Matlab-koden som implementerar DWO finns i appendix A.1
4.2
Lokal parameterskattning till direkt
viktopti-mering
Vi vill räkna ut prestandamåtten med (3.20), (3.21), (3.22) och (3.23) och sedan, som i (3.24), minimera med avseende på λ och använda det erhållna λ-värdet i DWO-metoden. Problemet med detta tillvägagångssätt är att det är svårt att mini-mera våra prestandamått med avseende på λ. Detta beror på att prestandamåtten använder information om vilka DWO-vikter som är noll, negativa eller positiva och det innebär att vi måste beräkna flera DWO-skattningar i förväg. DWO-metoden i sin tur innebär att ett QP som innehåller λ måste lösas. Vi kan alltså inte ändra λ tills vi hittar minimum utan att beräkna en eller flera DWO-skattningar för varje ändring.
Det vi kan göra är att välja några λ som vi vill undersöka och skapa våra prestandamått med dessa och sedan minimera med avseende på dessa valda λ-värden. Detta innebär dock att vi för det första kanske missar minimum eftersom vi inte tar med alla λ utan bara ett visst antal. För det andra måste ett stort antal DWO:er beräknas, vilket är mycket kostsamt. I fallet LCV måste vi dessutom räkna ut F1:N vilket innebär att ännu fler DWO:er måste beräknas.
Kapitel 5
Simuleringar
I detta kapitel kommer vi framförallt att diskutera ett antal simuleringar som har utförts för att undersöka prestandan hos de adaptiva metoderna vi får om vi använder prestandamåtten från avsnitt 3.2 på sättet som beskrevs i avsnitt 4.2 för att skatta en funktion f i punkterna ϕs(k), k = 1, . . . , Nϕs. Vi kommer att jämföra
dem med varandra och dessutom med vad DWO skulle ge för olika värden på λ i alla ϕs(k), k = 1, . . . , Nϕs. Några nya sätt att välja λ kommer också att införas.
För att skapa mätdata väljer vi först N stycken likformigt fördelade mätpunk-ter, ϕ(k), k = 1, . . . , N , mellan den minsta skattningspunkten ϕmin
s och den största
skattningspunkten ϕmax
s . För att få y(k), k = 1, . . . , N som i avsnitt 2.1 skapar vi
först f (ϕ(k)), k = 1, . . . , N , där f är en viss känd funktion och adderar därefter normalfördelat brus med standardavvikelse σ.
I tabell 5.1 anges några fler beteckningar som vi kommer använda för att be-skriva simuleringarna i detta kapitel.
5.1
Utan brus
Vi börjar med att testa DWO på data som saknar brus, alltså med σ = 0. Eftersom
σ = 0 måste vi använda metoden för σ = 0 som diskuterades i avsnitt 4.1.1, där
vi kan se att skattningen blir oberoende av L så länge L > 0 eftersom L bara kommer med i skalningen av målfunktionen. Att skattningen är oberoende av L innebär i sin tur att valet av λ inte är relevant.
5.1.1
Simulering 1: Utan brus
Vi börjar med att skatta en enkel sinus i ett simuleringsfall med värdena i tabell 5.2. Vi har valt L = 1 eftersom skattningen är oberoende av L. Resultatet visas i figur 5.1. Vi kan se att den skattade och den sanna kurvan är nästintill identiska.
Beskrivning
f (ϕ) funktionen som ska skattas
N antalet mätvärden
σ brusets standardavvikelse
λsann det sanna λ-värdet
λmax det största λ-värdet metoderna har att välja mellan
Nλ hur många λ-värden metoderna har att välja mellan
ϕmaxs den största skattningspunkten
ϕmin
s den minsta skattningspunkten
Nϕs antalet skattningspunkter
NI antalet iterationer i Monte Carlo-simuleringen
αmax det största α som används i LAIC och LFPE
ϕs(k), k = 1, . . . , Nϕs skattningspunkter
λval(k), k = 1, . . . , Nλ de λ-värden metoderna har att välja mellan
Tabell 5.1. Beskrivningar av beteckningarna som används i kapitel 5.
Värden Sim. 1 f (ϕ) sin(ϕ) N 100 σ 0 L 1 Nϕs 40 ϕs(k), k = 1, . . . , Nϕs 0.1667(k − 1)
Tabell 5.2. Funktion och värden som används vid simulering 1.
0 1 2 3 4 5 6 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 ϕ
Figur 5.1. Skattning med hjälp av DWO-metoden vid simulering 1. Skattningen
(hel-dragen), sanna funktionen (streckad) och mätvärden (stjärnor). Simuleringsvärden finns i tabell 5.2.
Värden Sim. 2 f (ϕ) sin(ϕ) N 100 σ 0.1 λsann 10 Nλ 21 Nϕs 40 ϕs(k), k = 1, . . . , Nϕs 0.1667(k − 1) λval(k), k = 1, . . . , Nλ 0.5(k − 1)
Tabell 5.3. Funktion och värden som används vid simulering 2.
0 1 2 3 4 5 6 −1 −0.5 0 0.5 1 ϕ
Figur 5.2. Skattningar med hjälp av DWO-metoden vid simulering 2. Skattningar för
de olika λval(k), k = 1, . . . , Nλ(heldragen), sanna funktionen (streckad) och mätvärden
(stjärnor). Simuleringsvärden finns i tabell 5.3.
5.2
Med brus
Nu när vi sett att DWO skattar funktioner bra utan brus kan vi börja undersöka hur DWO beter sig när det finns brus och låta metoderna från avsnitt 3.2 välja λ.
5.2.1
Simulering 2: Med brus och varierad inparameter
Innan vi låter de olika metoderna välja λ tittar vi på hur skattningen från DWO ändras med olika λ för att kunna välja våra λval(k), k = 1, . . . , Nλnågorlunda bra.
Vi fortsätter att skatta en sinus och i tabell 5.3 finns de valda värdena. I figur 5.2 ser vi hur skattningen ser ut för de olika λval(k), k = 1, . . . , Nλ. Den raka linjen
motsvarar noll och ju högre λ desto mer svängig blir skattningen. Vi ser att den största effekten av ett ändrat λ fås för låga λ-värden. Därför kommer vi från och med nu att välja våra λval(k), k = 1, . . . , Nλ logaritmiskt med start i 0.1, alltså
λval(k) = 10
−1+(k−1)(log10(λmax)+1)
0 2 4 6 8 10 12 14 16 18 20 0 1 2 3 4 5 6 7 8 9 10 k λ v a l ( k )
Figur 5.3. Logaritmiskt valda λval(k), k = 1, . . . , Nλmed Nλ= 20 och λmax= 10.
där λmax är det största av alla λval(k), k = 1, . . . , Nλ. Plotten i figur 5.3 visar ett
exempel på λval(k), k = 1, . . . , Nλ valt på detta sätt. Vi kommer också i
fortsätt-ningen välja ϕs(k), k = 1, . . . , Nϕs som Nϕs stycken ekvidistanta värden mellan
ϕmin
s och ϕmaxs .
5.2.2
Simulering 3: Låg andraderivata
I den tredje simuleringen har vi för första gången skattat λ med hjälp av LCV, LGCV, LAIC och LFPE. Eftersom två av våra metoder använder en designpara-meter som måste väljas i förväg kan ju detta liknas vid att välja ett λ direkt. Det kan alltså vara intressant att undersöka om det är lättare att välja α eller λ direkt. För DWO väljer vi ett konstant λ ∈ λvalför alla ϕs(k), k = 1, . . . , Nϕs.
Designpa-rametern α i LAIC och LFPE har också valts globalt för alla skattningspunkter och ligger i intervallet 2 ≤ α ≤ αmax för LAIC och 0 ≤ α ≤ αmax för LFPE.
Som mått på hur bra skattningen har varit väljer vi RMSE-felet från avsnitt 3.1.2, som vi kan använda tack vare att den sanna funktionen f är känd. RMSE-felet de olika metoderna ger kallar vi εLCV, εLGCV, εLAIC, εLF P E och εDW O. Vi
får här εLCV = s PNϕs k=1( ˆfLCV(ϕs(k)) − f (ϕs(k)))2 Nϕs ,
där ˆfLCV står för skattningen vi får med hjälp av LCV. Vi får εLGCV, εLAIC,
εLF P E och εDW O genom att byta ut ˆfLCV mot motsvarande funktionsskattning.
De valda värdena finns i tabell 5.4, i figur 5.4 ser vi vilka λ metoderna har valt och i figur 5.5 visas RMSE för olika val av designparametrarna i LAIC, LFPE respektive ren DWO. Som jämförelse visas också LCV och LGCV.
Vi ser i figur 5.4 att alla metoder förutom LFPE väljer λ mycket fluktuerande medan LFPE väljer mer konstant men nära det riktiga värdet. Vi ser också att LFPE verkar vara den enda metoden som inte hittar gropen i mitten. Det kan tyckas att detta borde göra att LFPE-skattningen blir sämst. Tittar vi på felet, som vi ser i figur 5.5, ser vi att det snarast är tvärtom. Det verkar vara bättre med
Värden Sim. 3 f (ϕ) sin(ϕ) N 100 σ 0.5 λsann 2 λmax 100 Nλ 20 ϕmin s 0 ϕmax s 6.5 Nϕs 40 αmax 50
Tabell 5.4. Funktion och värden som används vid simulering 3.
ett mer konstant värde eftersom LFPE är bäst av de vanliga metoderna. Notera också här att det α som ger det bästa slutresultatet för LFPE, någonstans runt
α = 18, är ett α som inte garanterar ett ickenegativt värde på måttet, se avsnitt
3.2.4. Vi kommer senare att se att i nästan alla fall blir det så att α > 2 ger det bästa slutresultatet.
I figuren kan vi också se att ren DWO med ett bra valt λ är ungefär lika bra som LFPE med ett bra valt α. LAIC i sin tur är bättre än både LCV och LGCV men bara för ett fåtal α. En annan observation är att LGCV är bättre än LCV trots sin enklare design.
5.2.3
Simulering 4: Hög andraderivata
En funktion med högre andraderivata än sinusfunktionen från de tre första simule-ringarna har också undersökts och de tillhörande värdena finns i tabell 5.5. Detta är samma funktion och värden som i exempel 7.1 i Roll (2003) men med annan slumpmässig data. Om vi tittar i figur 5.6 ser vi att observationerna om valda λ från förra simuleringen kan göras här också även om nedgången i mitten ligger något snett. Om vi tittar på RMSE i figur 5.7 ser vi att nu är LCV och LGCV i det närmaste exakt lika bra medan ren DWO är ganska klart bättre än LFPE och LAIC aldrig blir bättre än LCV och LGCV.
5.3
Monte Carlo-simuleringar
Att titta på enskilda simuleringar ger inte så mycket i vårt fall eftersom relativt mycket beror på slumpmässiga data. Det är därför naturligt att använda återupp-repning av likadana simuleringar för att lättare kunna dra slutsatser från dem. Detta brukar kallas Monte Carlo-simuleringar.
0 1 2 3 4 5 6 10−1 100 101 102 ϕ (a) LCV 0 1 2 3 4 5 6 10−1 100 101 102 ϕ (b) LGCV 0 1 2 3 4 5 6 10−1 100 101 102 ϕ (c) LAIC (α = 4) 0 1 2 3 4 5 6 10−1 100 101 102 ϕ (d) LFPE (α = 18)
Figur 5.4. Valda λ (heldragen) för olika metoder och |f00(ϕ)|/σ (streckad) vid simulering
3. Simuleringsvärden finns i tabell 5.4
0 20 40 60 80 100 0.2 0.25 0.3 0.35 0.40 10 20 30 40 50 0.2 0.25 0.3 0.35 0.4
Figur 5.5. RMSE-felen för olika metoder vid simulering 3. εLCV (heldragen horisontell),
εLGCV (streckad horisontell), εLAIC(streckad), εLF P E(heldragen) och εDW O(heldragen
med prickar). Simuleringsvärden finns i tabell 5.4. Övre x-axeln anger vilket värde på α som använts i LAIC och LFPE. Undre x-axeln anger vilket värde på λ som använts i den rena DWO:n.
Värden Sim. 4 f (ϕ) ϕ2sin(2ϕ) N 100 σ 1 λsann 12.87 λmax 20 Nλ 20 ϕmin s -2 ϕmax s 2 Nϕs 40 αmax 50
Tabell 5.5. Funktion och värden som används vid simulering 4.
−2 −1 0 1 2 100 101 ϕ (a) LCV −2 −1 0 1 2 100 101 ϕ (b) LGCV −2 −1 0 1 2 100 101 ϕ (c) LAIC (α = 2) −2 −1 0 1 2 100 101 ϕ (d) LFPE (α = 23)
Figur 5.6. Valda λ (heldragen) för olika metoder och |f00(ϕ)|/σ (streckad) vid simulering
0 5 10 15 20 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 10 20 30 40 50 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figur 5.7. RMSE-felen för olika metoder vid simulering 4. εLCV (heldragen horisontell),
εLGCV (streckad horisontell), εLAIC(streckad), εLF P E(heldragen) och εDW O(heldragen
med prickar). Simuleringsvärden finns i tabell 5.5. Övre x-axeln anger vilket värde på α som använts i LAIC och LFPE. Undre x-axeln anger vilket värde på λ som använts i den rena DWO:n.
5.3.1
Simulering 5, 6 och 7: Låg andraderivata och
varie-rande brusnivå
Vi börjar med en funktion som har lågt L-värde skattad med varierande brusni-vå. Värdena finns i tabell 5.6, där NI står för antalet iterationer i Monte
Carlo-simuleringen.
Vi inför ¯εLCV, ¯εLGCV, ¯εLAIC, ¯εLF P Eoch ¯εDW Osom är medelvärdena av εLCV,
εLGCV, εLAIC, εLF P E och εDW O över de olika iterationerna. I figur 5.8 ser vi att
om vi inte vill välja något värde i förväg är LCV och LGCV ungefär lika bra men att LFPE och även ren DWO är mycket bättre för många val av respektive designparameter. I simulering 6 är till exempel LFPE bättre än LCV och LGCV för alla α mellan ungefär 10 och 40 och den rena DWO:n är bättre än LFPE:s bästa värde för λ mellan ungefär 0.5 och 1. Den enda gången LAIC verkar värd att använda är i simulering 7 då den ger ett bra resultat för de flesta α. Det verkar vara så att LAIC blir bättre relativt de andra metoderna för ökat brus.
Vi vet sedan tidigare att valet av λ kan ses som ett mått på hur snabbt vi tillåter skattningen att variera. Detta innebär att om vi vill dämpa ett stort brus ska ett lägre λ väljas men om ett för lågt λ väljs tillåter vi inte skattningen att variera tillräckligt snabbt för att kunna följa med den sanna funktionen. Storleken på λ blir alltså en avvägning mellan felet som ges av att modellen inte kan återskapa den sanna funktionen även om bruset skulle vara noll (s.k. biasfel, se t.ex. Ljung och Glad (2004)), i vårt fall för lågt λ så att skattningen inte tillåts svänga tillräckligt mycket, och felet som bruset ger (s.k. variansfel, se t.ex. Ljung och Glad (2004)). Även α i metoderna LAIC och LFPE har denna egenskap men i motsatt riktning, alltså ett för litet α ger variansfel och ett för stort ger biasfel.
Om vi tittar i figur 5.8 igen ser vi att det blir som väntat, man ska välja lägre
Värden Sim. 5 Värden Sim. 6 Värden Sim. 7
f (ϕ) sin(ϕ) sin(ϕ) sin(ϕ)
N 100 100 100 σ 0.5 1 2 λsann 2 1 0.5 λmax 10 10 10 Nλ 20 20 20 ϕmin s 0 0 0 ϕmax s 4π 4π 4π Nϕs 100 100 100 αmax 50 50 50 NI 50 50 50
Tabell 5.6. Funktion och värden som används vid simulering 5, 6 och 7.
0 2 4 6 8 10 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50 10 20 30 40 50 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 (a) Simulering 5 0 2 4 6 8 10 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.50 10 20 30 40 50 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 (b) Simulering 6 0 2 4 6 8 10 0.5 0.6 0.7 0.8 0.9 10 10 20 30 40 50 0.5 0.6 0.7 0.8 0.9 1 (c) Simulering 7
Figur 5.8. RMSE-felen för olika metoder vid simulering 5, 6 och 7. ¯εLCV (heldragen
ho-risontell), ¯εLGCV (streckad horisontell), ¯εLAIC (streckad), ¯εLF P E(heldragen) och ¯εDW O
(heldragen med prickar). Simuleringsvärden finns i tabell 5.6. Övre x-axeln anger vilket värde på α som använts i LAIC och LFPE. Undre x-axeln anger vilket värde på λ som använts i den rena DWO:n.
