Dag 3

Report
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-163.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).

similar documents