Transcript Dag 3
Hvordan får man data og modell til å passe sammen? Ekstremverdi-analyse Målet er å estimere T-års-ekstremen (flommen). T-års-ekstremen er slik at etter T år vil det i forventning være én overskridelse av T-års-ekstremen. For årlige ekstremer blir dette 1/T-kvantilen til fordelingen disse verdiene. Data: Enten maksima/minima fra blokker eller fra maksima/minima over/under terksel. For maksima/minima fra årsblokker blir dette klassisk gjort ved å estimere ekstremverdi-fordelings-parameterne, og hente 1/T-kvantilen derifra. Asymptotisk teori angir standardfordelinger hvis man har et sett maksima over gitte (store) blokker med uavhengigeFordeling (svart), data (søyler), MLdata (GEV) eller over en gitt stor verdi estimert fordeling (rød), Bayesiansk (Pareto). prediksjonsfordeling (blå). Ekstremverdi-analyse-problemer Merk at sannsynligheten for å overstige en 100-års-flom på en 100- års-periode ikke er 100%. Hvert år en sannsynligheten for overstigning 1/T, som over T uavhengige år blir ca. lik 1-e-163.2%. Merk at forutsetningene for asymptotikken er brutt i NVE-data (slettes ingen uavhengighet innenfor år). Merk at estimert T-års-ekstrem vil være det vi får fra å velge en fordelingsfamilie, estimere parametre og beregne 1/T-kvantil fra dette. Dette er ikke den egentlige T-års-ekstremen, fordi vi er usikre på korrekt fordelingsfamilie og usikre på parameterverdiene. Dette kan føre til skjevheter i estimert T-årsekstrem. Eks: Trekker man et datasett på 25 år trukket fra en Gumbel-fordeling med tilfeldige parametre, vil ML-estimert 100-års-flom overstiges en av 65 ganger når man trekker en ny årsmaks fra Gumbel-fordelingen. I praksis oppfører derfor estimert 100-års-flom seg som en 65-års-flom. Det samme fås for L-momenter, men ikke Bayesianske prediksjonsfordelinger tatt fra vag prior. Regresjon Regresjon er når en stokastisk variabel (respons) antas å avhenge av andre variable (forklaringsvariable, som i denne sammenheng ikke antas være stokastiske). En del av variasjonen en ser i respons-variabelen er altså forklart via variasjon i andre variable. vekt Eksempel: Vekt (respons) versus høyde (forklaringsvariabel) høyde Lineær regresjon En lineær regresjon, undersøker vi en lineær sammenheng mellom respons og forklaringsvariable: Y=0+1x1+2x2+…+pxp • Merk at modellen er lineær i koeffisientene, 0,…,p, ikke nødvendigvis i forklaringsvariablene. Så modellen Y= 0+1x+2x2 er en lineær modell. • Den statistiske modellen bak er som følger: Yi 0 1 x1, i 2 x 2 , i ... p x p , i i er uavhengig støy. der i ~ N ( , ) Lineær regresjon - eksempel vekt Eksempel: vekt i a b * høyde i i Regresjonskoeffisientene, a og b, kan tilpasses via ML-estimering. Grafen viser en slik tilpasning. Regresjonen ser ut til å beskrive det som er å skimte av systematikk i dataene. Modellen selv er dog snål. Ifølge den skal det finnes en høyde slik at du kan forvente null vekt samt at du via tilfeldigheter kan ha negativ vekt selv der en forventer positiv vekt (dette pga normalfordelings-antagelsen). Man kan redde denne situasjonen ved å anta at det er log-transformert høyde og vekt som kan beskrives via lineær regresjon. Dette betyr en power-law for originaldata. vekt i A * høyde b høyde vekt der E i ~ log N ( 0 , ) i * Ei høyde Lineære regresjon – når man går amok i forklaringsvariable Med de muligheter som ligger i regresjon, kan man falle for fristelsen til å bare legge på flere og flere forklaringsvariable. Som et eksempel, kan vi legge på vekt høyere-ordens polynomledd i høydemot-vekt-eksempelet: 2 3 4 v i 0 1 hi 2 hi 2 hi 4 hi i Det som skjer er at tilpasningen til data blir bedre (alltid!), men en kan forvente at evnen til å forutse utkommet av nye data (prediksjon) blir bare verre. Sammenhengen selv blir mer kaotisk og parameter-usikkerhetene blir større og større. Dermed blir prediksjonsusikkerheten større. høyde Hvordan unngå å gå amok? Det er i basis to muligheter for å unngå å gå amok i forklaringsvariable. 1. Tenk gjennom dataenes natur (som betyr power-law heller enn lineærmodell for vår vekt-mot-høyde) og hva du ønsker å gjøre med din regresjon. 2. Bruk hypotesetesting (modellvalg) til å begrense deg. (PS: R rapporterer p-verdier for alle forklaringsvariable). Det siste kan gjøres ved å: a) b) c) d) Starte med en enkel modell og legge til variable så lenge du finner noen som er statistisk signifikante Starte med en tilstrekkelig komplisert modell og ta vekk variable så lenge de ikke er signifikante. Gå igjennom alle tenkelige modeller og velg den med best informasjonskriterie. (Ikke anbefalt for store antall forklaringsvariable!) Bruke Bayesiansk metodikk. Merk at i høyde-vs-vekt-eksempelet er ikke høyde signifikant i utgangspunktet! Usikkerhet Estimatorene i regresjon kommer med en viss usikkerhet. Disse blir rapportert i R. vekt Prediksjons-usikkerhet Estimasjonsusikkerhet Når konfidensintervallene omslutter 0, betyr det at en ikke kan forkaste at en forklaringsvariabel har null effekt. M.a.o. at den er ikke statistisk signifikant. Dette påvirker usikkerheten estimatet for den virkelige sammenhengen mellom respons og forklaringsvariable, altså forskjellen mellom Y 0 1 x1 2 x 2 ... p x p og Yˆ ˆ 0 ˆ1 x1 ˆ 2 x 2 ... ˆ p x p samt usikkerheten til prediksjoner: Yˆ pred ˆ 0 ˆ1 x1 ˆ 2 x 2 ... ˆ p x p Prediksjoner er mer usikre enn estimat, siden man i tillegg får de individuelle variasjonene på toppen av estimasjons-usikkerhetene. høyde Simulert datasett Residualer Residualer er avviket mellom måling og modell på y-aksen (responsen). Disse avvikene kan si noe om hvorvidt modellantagelsene er riktig. 1. 2. 3. En tydelig trend i residualene antyder at funksjonssammenhengen kan være gal. Er trenden i tid, tyder det på at gradvis forandring i umålte forklaringsvariable spiller en rolle eller at man har å gjøre med korrelasjon i tid (tidsserie). Hvis residualene ikke later til å normalfordelt, kan det tenkes en transformasjon trenges, eller at en annen type regresjon er nødvendig. Også hvis variasjonen i residualene har en trend (”trumpetform”), er støyleddene modellert feil (heteroskedastisitet). Remodellering (mer avansert regresjon) eller datatransformasjon kan være nødvendig. Data+ regresjon Data+ regresjon residualer residualer qq-plott Data+ regresjon residualer Ikke-lineær regresjon Ikke all regresjon er lineær. Noen ganger trenger vi å lete etter sammenhenger mellom respons og forklaringsvariable som har en annen form. Et eksempel er vannføringskurve-tilpasning med ukjent bunnvannstand: Q=C(h-h0)b Selve etter en log-transformasjon, ødelegger h0 lineariteten: q=a+b*log(h-h0) ML-optimering er fremdeles mulig, men kun via numeriske metoder. F.eks. i vf-kurve-tilpasning vil man kunne optimere parametrene a og b analytisk, men h0 må optimeres numerisk. For mer kompliserte modeller, kan sofistiske optimeringsmetoder bli nødvendig. (Evt. Bayesianske metoder.) En fare med kompliserte ikkelineær modeller er at likelihood’en kan ha flere topper (multimodalitet). Vannføringskurvetilpasning på Gryta Skal nå se på Gryta stasjon, uten å anta at h0=0. Vi vil bruke ”brute force” ved å se på et intervall av mulige h0-verdier fra minste målte vannstand, hm, til hm-100m. Ser ut som vi kan maksimere loglikelihood (og dermed også likelihood) med en verdi for h0 nærme null. En nærmere titt gir optimal h0=+8cm. Merk de tidligere nevnte urimelige parameter-estimatene som av og til kan oppstå. Bayesiansk regresjon Skal igjen se på Gryta stasjon. Under Bayesiansk regresjon antas en førkunnskap. Denne kan trekkes fra samlingen av norske stasjoner, men for stasjonen Gryta vet vi at nullvannstanden ligger rundt h0=0 og siden det er et V-overløp vet vi også at b ca. lik 2.5 bør være en grei hydraulisk antagelse. I VFKURVE3 settes førkunnskapen i et eget vindu. Merk at Bayesiansk statistikk har mindre problemer med å håndtere multimodalitet. Simulering fra a’ posteriori-fordelingen blir dog vanskeligere, men det finnes dog relativt effektive metoder for å håndtere dette. Bayesiansk regresjon 2 Man foretar så analysen, som vil trekke masse parametre fra a’ posteriorifordelingen. I tillegg til å gi estimater, gir dette også en pekepinn på parameterusikkerheten. For parametre der vi satt en skarp førkunnskap, vil typisk a’ posteriorirfordelingen være innenfor det skarpe intervallet. Siden vi får oversikt over parameterusikkerheten vil også kurve-usikkerheten være tilgjengelig på fordelings-form. Med mye data og/eller bra førkunnskap, kan kurveusikkerheten bli svært liten. Regresjon mellom tidsserier Hvis vi ønsker å kjøre regresjon av en vannføringsserie mot en annen, havner vi på litt dypt vann, siden modellantagelsene ikke er tilstede (avhengighet i støyleddene). Teorien sier dog at estimatene vil være forventningsrette. Men usikkerhet og modelltesting vil bygge på antagelser som kan være radikalt feile. Typisk vil usikkerheten bli sterkt undervurdert. Her er to uavhengig simulerte tidsserier. Plotter vi den ene mot den andre, kan det se ut som det er en hvis avhengighet, noe en lineær regresjon vil støtte. Men dette skyldes kun at begge seriene har tidsavhengighet! Resultat fra R, summary(lm(x2~x1)): x1 -0.47232 0.04747 -9.95 < 2e-16 *** Tidsserie-analyse Statistiske tidsserier er data i tid, der det en eller annen form for avhengighet mellom det som skjer på et tidspunkt og det som skjer i neste. Eksempel: vannføringsserier, magasinering, nedbør for fin tidsoppløsning… Hvis tidsavhengigheten ikke tas hensyn til, vil man svært ofte undervurdere usikkerhetene involvert og man kan ikke stole på utfallet av modelltesting. Når modell krasjer med virkelighet 3 – uavhengig støy vs tidsserie Har simulert “vanntemperatur” med forventing =10. Antar kjent varians, =2. Ønsker å estimere og teste =10. Modell 1, avhengighet: Ti=+i, i~N(0,1) u.i.f. - Grafen ser ut til å fortelle en annen historie... - Estimert: ˆ x 11 . 4 , sd ( ˆ ) s / n 0 . 2 - 95% konf. int. for : (11.02,11.80). =10 forkastet med 95% konfidens! • Modell 2, auto-korrelert modell med forventning , standardavvik og auto-korrelasjon a. – Lineær avhengighet mellom temperaturen en dag og neste. – Estimert: ˆ x 11 . 4 , sd ( ˆ ) 1 . 4 – 95% konf. int. for : (8.7,14.10). =10 ikke forkastet. Tidsserier – diagnostiske plott Det er flere måter å få innblikk i en tidsseries natur. 1. Autokorrelasjon. Dette er et plott som viser korrelasjonen mellom verdien på et tidspunkt og et gitt antall tidskritt videre, som funksjon av disse tidssskrittene. Normalt vil dette avta etter hvert, men for serier med sesongavhengighet, kan det hende du får en negativ avhengighet etter et halvår og en ny positiv avhengighet etter et helt år. 2. Fourier-analyse. Dette dekomponerer en tidsserie inn i sinus/cosinusfunksjoner med ulik periodisitet. Tidsserier med sesong-avhengighet vil da ha en sterk topp på ett år. Diagnostikk og sesong-avhengighet For mange hydrologiske tidsserier vil sesong-avhengighet være opplagt. Men hva er tidsserienes natur etter at man har tatt hensyn til dette? I start-systemet er det en opsjon kalt ”konform transformasjon” som trekker fra årsgjennomsnittet og deler på standardavviket. Dermed kan autokorrelasjon ses når sesongavhengigheten er (mer eller mindre) tatt vekk. Uten en slik operasjon, vil en analyse på temperaturdata typisk angi en korrelasjonstid (tid før korrelasjonen går under en viss grense, som for eksempel 0.5) på opptil flere år. Etter operasjonen, vil en typisk korrelasjonstid være på rundt en uke. Altså, hvis man tar hensyn til sesongenes svinginger, er dagens temperatur kun en pekepinn på fremtidens temperatur rundt en uke frem i tid. Statistiske tidsseriemodeller Det finnes et arsenal av statistiske tidsserie-modeller. En stor gruppe av disse, kalles ARIMA modeller. Dette er sammensatt av kombinasjoner av modeller som har følgende elementer: AR (autoregressive) I (integrerte) og MA (moving average). AR-modeller: Dette er modeller der neste verdi avhenger av en gitt mengde av de foregående verdiene. F.eks. AR(1) avhenger kun av siste verdi, som er det som er kjent som en Markov-kjede: X t X t 1 (1 ) t der t ~ N ( 0 ,1) er uavhengig støy MA-modeller: Modeller basert på glidende midling: X t t 1 t 1 ... p t p der t ~ N ( 0 , ) er uavhengig støy Integrerte modeller: Dette er modeller der man transformerer data fra originaltidsserien til differanser: Yt X t X t 1 Dette gjøres for å modellere tidsserier som ikke er stasjonære, dvs. som ikke har noe fast fordeling eller forventningsverdi. Mer diagnostikk En MA-modell vil gi autokorrelasjonsplott (acf) som brått dør hen. Dør den hen etter k tidskritt, har man å gjøre med en MA-k-modell. En AR-modeller kan undersøkes ved et tilsvarende plott kalt ”partial autocorrelation function” (pacf). Data produsert av en AR-k-modell vil ha et pacf plott med bare k signifikante verdier i starten. Her et eksempel på et pacf-plott, tatt på data generert fra en AR(1)modell: Kalman-filter Et Kalman-filter er basert på en modell som har en skjult tidsseriene styrt av en multidimensjonal AR(1)-prosess. På toppen av disse har man observasjonene. Merk at dette rammeverket kan brukes til å binde sammen flere tidsserier i en ”tidsserie-regresjon”. tid Tilstand: Observasjoner: X1 X2 X3 Xn Y1 Y2 Y3 Yn For lineære modeller er dette en metode som analytisk er i stand til å regne ut forventing og varians for de skjulte tidsseriene betinget på observasjonene, samt for normalfordelte variable å regne ut likelihood. En modell med skjulte tilstander som skal infereres mhp observasjoner, er i stand til å håndtere manglende data. Dette kan dermed passe bra til utfylling av hull i tidsserier. Eksempel på bruk av Kalman-filter I dette eksemplet blir tre temperaturserier nær hverandre brukt. En del data ble fjernet og et Kalman-filter med korrelert støy-ledd mellom de tre seriene, ble undersøkt. Plottene viser ifyllingen av manglende data, samt usikkerhet og de dataene som ble tatt vekk. Siden modellen tillater korrelasjoner, vil data fra en stasjon informere om hva som skjer en annen plass. Der det mangler data på flere stasjoner, vil usikkerheten ”boble” ut. Kontinuerlig-tid stokastiske prosesser Selv om målinger gjøres på diskrete tidspunkt, er det vi henter data fra gjerne kontinuerlig i tid (vannføring, f.eks). Ofte kan tidsoppløsningen forandre seg også, som gjør at modeller med diskret tid feiler. Kont. tidsserie-modeller gir en sannsynlighet for fremtidige utfall på vilkårlige tidspunkt, gitt historikken. Kan også brukes til interpolasjon. Gir usikkerhet så vel som estimat. Eksempler: Poisson-prosessen (diskrete hendelser i kontinuerlig tid) Birth-death-prosesser (antallsdata i kontinuerlig tid) Wiener-prosessen (random walk) Ornstein-Uhlenbeck Lagdelte lineære modeller Stokastiske diffligninger Levy-prosesser t t t 1.96 s t -1.96 s Skjult OU Målt prosess påvirket av OU Romlige modeller og tid-roms-felt Interpolasjon og ekstrapolasjon er noe som kan være aktuelt i rom så vel som i tid. Har man et forhold til romlige avhengigheter, kan man bruke statistikk til å gjøre slik type estimering og si noe om estimeringsfeilen. Modellene kan være diskret eller kontinuerlig. Ofte brukt metodikk, ”kriging”, som antar en funksjonsform på avhengighetsstrukturen (semi-variogram) og kjører en regresjon på denne funksjonen mot estimerte avhengighetsmål. Alternativ: ML eller Bayesiansk analyse på avhengighetsstruktur. INLA. Utvidelse: Tid-roms-felt, altså avhengighetsstrukturer i både tid og rom. For å fylle ut en funksjon i både tid og rom (som f.eks. nedbør eller temperatur).