Transcript Powerpoint

Hvordan får man data og modell til å passe sammen?
Faget statistikk
 Ikke tørre tall, men essensen i dem.
 Modell vs data – tilpasning av interesse-parametre
 Eks på parametre: gjennomsnittelig årsnedbør, en
vannføringsseries forklaringsverdi på en annen slik serie,
magasinering som funksjon av nedbørsareal.
 Parametre er i utgangspunktet ukjent, men dataene gi oss
et estimat samt en antydning om hvor usikre disse
estimatene er.
 Modellvalg – gir svar på spørsmål
 Eks: Er årsnedbøren lik i to nabofelt? Kan vi si noe som helt
om en vannføringsserie på bakgrunn i en annen?
 Svarene er ikke absolutte, men gis med en viss sikkerhet.
Datausikkerhet
 Perfekte målinger + perfekte modeller =
Null datausikkerhet =>
Ingen parameter- eller modell-usikkerhet
Datausikkerhet:
 Reelle målinger er dog beheftet med usikkerhet.
 Modellene kan ikke ta alt med i betraktningen. Umålte
”confounders” (lokal topografi og jordsmonn i en
hydrologisk modell, for eksempel.)
 Håndteres ved måten målingene sprer seg på, i.e.
sannsynlighetsfordelingen.
Modellusikkerhet
 Datausikkerhet => Usikkerhet i parameterverdier og
hvilken modell som er best (modellusikkerhet)
Håndtering av modellusikkerhet:
I. Usikkerhet i parameterverdier og modellvalg kan
håndteres likt som data: ved
sannsynlighetsfordelinger. Bayesiansk statistikk.
II. Parameterverdier og modeller kan håndteres som
fikserte men ukjente – frekventistisk statistikk
(klassisk).
Sannsynlighet
 Sannsynlighet:
I.
II.
III.
Angir langtidsraten av utfall som havner i en gitt kategori. F.eks. vil
1/6 av alle terningkast gi utfallet ”en”.
Angir forholdet mellom en gevinst og hva du er villig til å risikere for
den. F.eks. kan du være villig til å risikere 10kr for å få tilbake 60kr
hvis du får ”en” på en terningkast.
Kan gi en formell beregningssystem for usikkerhet og forventning.
Sannsynlighet 1/6 for å få ”en” på et terningkast antyder at du ikke
har noen større eller mindre grunn til å forvente ”en” enn noe annet
utfall på terningen.
II og III er begge Bayesianske sannsynligheter, som kan oppfattes som
“subjektive” mens I erfrekventistisk og “objektiv” i den forstand at sannsynlighetene
antas komme fra tings iboende egenskaper. (Så spørs det om dette virkelig er
tilfelle.)
Sannsynlighetlover
1. For en hendelse A skriver vi
sannsynligheten for
hendelsen som Pr(A).
2. 0≤Pr(A)≤1
3. Pr(A)+Pr(ikke A)=1
Eks:
Pr(”Du får en ener på ett
terningkast”)
Pr(flom på vestlandet)=1.1
betyr at du har regnet feil.
Pr(”to eller mer på et
terningkast) = 1-Pr(”ener”)
= 1-1/6=5/6
4. Pr(A eller B)=Pr(A)+Pr(B) når Pr(”ener eller toer”) =
Pr(”ener”)+
A og B ikke kan stemme
Pr(”toer”)=1/6+1/6=1/3
samtidig.
Sannsynlighetlover 2 – betinget
sannsynlighet
Eks:
Pr(A | B) gir sannsynligheten for A
under forutsetning at B stemmer.
Pr(A|B)=Pr(A) betyr at A og B er
avhengige. B gir ikke informasjon
om A. Motsatt gir B informasjon om A,
Hvis andre terningkast ikke lar seg
påvirke av første, vil
Pr(”ener på andre” | ”ener i første”) =
Pr(”ener på andre”).
som er drivkraften i Bayesiansk statistikk. Lar vi B=”ener i første kast” og
A=”ener i første kast”:
Pr(A og B)=Pr(A | B)Pr(B) der Pr(A|B) Pr(”ener på første og andre terningkast”)
= Pr(A|B)Pr(A) = Pr(A)Pr(B) =
betyr ”sannsynligheten for A gitt B”. 1/6*1/6=1/36.
Siden Pr(A og B)=Pr(B|A)Pr(A) også, Fra Bayes teorem: Hvis B ikke gir
informasjon om A, Pr(A|B)=Pr(A), så
får vi Bayes formel:
gir ikke A informasjon om B heller,
Pr(A|B)=Pr(B|A)Pr(A)/Pr(B)
Pr(B|A)=Pr(B).
Sannsynlighetsfordelinger –
endelige utfall
En sannsynlighetsfordeling gir et hvert mulig utfall en
sannsynlighet.
Eks: En terning
Alle utfall fra en til seks er like
sannsynlige
Sum av to terninger
En sum på tre (2+1 eller 1+2) er dobbelt så
sannsynlig som et utfall på 2 (1+1).
Sannsynlighetsfordelinger –
kontinuerlige utfall
f(x)=sannsynlighetsfordeling
En sannsynlighetsfordeling med kontinuerlige utfall gir
et hvert mulig intervall en sannsynlighet. Dette heter
gjerne en sannsynlighetstetthet.
Eks: uniform fordeling:
I dette tilfelle er f(x)=1 for alle x mellom 0 og en og
0 utenfor.
Hva dette sier, er at utfall mindre enn 0 eller større
enn 1 er umulig. Videre sier det at alle intervaller
innenfor (0,1) som har lik størrelse, er like
sannsynlige. Sannsynligheter må summeres til en
og sannsynligheten for to ulike utfall er summen
av enkeltsannsynlighetene. Dermed blir
sannsynligheten for et utfall i et intervall
proporsjonalt med størrelsen til intervallet.
1
0
1
x=utfall
Normalfordelingen
Til forskjell fra uniform fordeling er alle utfall på tallinjen mulig, men den
har likevel et klart senter og en klar utspredning.
Toppen angir både fordelingens ”gjennomsnittsverdi”,
forventningsverdien, samt stedet der det er 50% sannsynlighet for at
utfallene havner over og 50% for at de havner under (medianen).
Spredningen, kjent som standardavviket,
angir et område der det er mest
sannsynlig å finne utfallene. Det er 68%
sannsynlighet for at utfallet er ett
standardavvik unna forventningsverdien.
Det er 96% sannsynlighet for at utfallet
er innenfor to standardavvik.
Sannsynligheten avtar raskt utover.
Forventning og standardavvik er to parametre som til sammen
spesifiserer normalfordelingen.
Mer om normalfordelingen
Fordelingsfunksjonen, f(x), er glatt. Sannsynligheten for å få et utfall i et
lite intervall (x,x+dx) er f(x)*dx.
Matematisk ser den slik ut:
der  er forventingsverdien og  er
standardavviket.
 ( x   )2 
1

f ( x) 
exp 
2
2 
2 

Skal man regne ut sannsynligheten for å
få et utfall i et vilkårlig stort intervall må
man summere sannsynligheten for
masse små. En slik sum er kjent som et
integral.
At en tilfeldig (stokastisk) variabel, X, er
normalfordelt, skriver vi som: X~N(,).
Hvorfor normalfordelingen?
Selv om normalfordelingen ser litt komplisert ut matematisk, har den en
rekke gode egenskaper.
 Den er glatt og tillater alle mulige utfall.
 Er karakterisert med en enkelt topp.
• Det viser seg at hvis du betinger på at en funksjon er positiv, glatt og har
bare en topp, vil normalfordelingen være den enkleste og en som lokalt
tilnærmelsesvis er lik enhver annen fordeling med samme egenskaper.
 Symmetrisk
 Informasjonsmessig er det den fordelingen som koder for en gitt
sentrering (forventning) og spredning (standardavvik) med minst
mulig ekstra informasjon. (Maksimal entropi).
 Summen av to normalfordelte størrelser er normalfordelt.
 En stor sum av størrelser med lik fordeling vil være ca.
normalfordelt. (Sentralgrenseteoremet).
 Matematisk behagelig å jobbe med (tro det eller ei!)
Skalastørrelser – lognormalfordelingen
Når en størrelse er nødt til å være strengt positiv (massen til en person,
volum i et magasin, vannføringen i en elv), passer det ikke å bruke
normalfordelingen.
En enkel måte å fikse dette på, er å ta en logaritmisk transformasjon på
størrelsen. Hvis en stokastisk variabel X>0, vil log(X) anta verdier
over hele tall-linjen.
Antagelsen log(X)~N(,) gir også en fordeling for X, kalt den
lognormale fordelingen, X~log N(,).
 (log(x)   ) 2 
1

f ( x) 
exp 
2
2
2 x


Ekstremverdifordelinger
Ekstremverdifordelinger er fordelingstyper som typisk vil være gode
tilnærmelser til fordelingen til ekstreme hendelser, under gitte
betingelser. Betingelsene vil angi hvilken fordeling det er snakk om.
1. Maksimum/minimum over et
gitt tidsintervall. Eks: årsflommer
Her sier teorien det er GEVfordelingen som gjelder. Denne har
tre parametre, en som angir
sentrering, en for spredning og en
angir formen.
f ( x) 
1

t ( x) 1 e t ( x )
x   1/ 


))
når   0 
(1   (
der t ( x)  


( x   ) / 


e
ellers


Ekstremverdifordelinger 2
1. Maksimum over en gitt terskelverdi
Her sier teorien det er Pareto-fordelingen som gjelder. Denne har to parametre, en
som angir nedre grense, xm, og en som angir formen, .
f ( x) 
 xm 
x
 1
for x  xm
Pareto-fordelingen kan være
ekstremt tunghalet, det vil si at
sannsynlighets-tettheten avtar
veldig lite utover. Dette kan være
problematisk for forventing og
standard-avvik.
Egenskaper til stokastiske variable kvantiler
 Ut ifra fordelingsfunksjonen, kan man få sannsynligheten for å ligge
inne i et intervall eller over en gitt verdi eller under en gitt verdi.
 Snur man på dette, kan man få en kvantil/persentil, q(p). Dette er en
verdi slik at sannsynligheten for at X skal ligge under denne er p.
 En spesialversjon av dette er medianen, der det er 50% sannsynlighet
for at X skal anta verdier over median-verdien og 50% for at X skal anta
verdier under. Eks: Hvis X~N(,), så er medianen .
 Kvantiler kan brukes til å angi hva som kan antas være rimelige utfall.
Ofte brukes 2.5%- og 97.5%-kvantilen. Det er 95% sannsynlighet for at
X skal være innenfor dette intervallet. Eks: hvis X~N(,), så er X
innenfor (-1.96 , +1.96) med 95% sannsynlighet.
Egenskaper til stokastiske variable forventingsverdi
 Forventningen er en stokastisk variabels
”gjennomsnittsverdi”. For et endelig utfall, finner man dette
ved å vekte et utfall med dens sannsynlighet og summere
over mulige utfall.
N
E ( X )   xi Pr(X  xi ) der det er N ulike mulige ut fall,(x1 ,..., xN ).
i 1
 Eks:
i.
ii.
iii.
iv.
v.
For en terning er forventningsverdien 3.5.
For en uniformt fordelt variabel mellom 0 og 1, er forventingen ½.
For en normalfordelt variabel er forventingen .
For en lognormalfordelt variabel er forventingen exp(+2/2)
En Pareto-fordelt variabel har ikke forventing for <1.
Egenskaper til stokastiske variable
– standardavvik og varians
 Standardavviket angir hvor mye en stokastisk
variabel sprer seg på. Teknisk sett er den
kvadratroten av variansen, som er forventet
kvadratisk avvik fra forventingsverdien:
N
Var( X )   ( xi  E ( x))2 Pr(X  xi ) der det er N ulike mulige utfall,(x1 ,..., xN ).
i 1
 Eks:
ii.
For en uniformt fordelt variabel mellom 0 og 1, er variansen 1/12.
For en normalfordelt variabel, er standardavviket .
iii.
En Pareto-fordelt variabel har ikke varians eller standardavvik for <2.
i.
Trekninger av stokastiske variable
 Hvis vi er i stand til å trekke fra en statistisk
fordeling, vil vi med mange nok trekninger se at:
i.
Gjennomsnittet nærmer seg forventingsverdien.
1 n
x   xi  E ( X )
n i 1
ii.
Kvadratavviket
nærmer seg variansen.
n
S2 
iii.
1
( xi  x ) 2  Var( X )

n  1 i 1
Histogrammet nærmer
seg fordelingsfunksjonen
(trekningenes kvantiler nærmer
seg fordelingens kvantiler.)
f(x)
Diagnostikk på
fordelingsfunksjoner
 Man kan vise histogrammet til
dataene og sammenligne med
fordelingen.
 Eventuelt kan man plotte
teoretiske kvantiler mot datakvantiler, såkalte qq-plott. Har
man rett fordeling, skal disse
kvantilene ligge på en rett linje.
Statistisk inferens
 I realiteten kan det være usikkerhet om hvilken fordeling
(modell) som passer til å beskrive hvordan dataene har blitt
produsert.
 Gitt modellen, vil likevel parameterverdiene være ukjent.
Naturen vil ikke bare dumpe dette i hendene våre.
 Statistisk inferens dreier seg om å bruke data til å si noe om:
Modellvalg
ii. Usikkerhet rundt modellvalget
iii. Estimering av parameterverdier
iv. Usikkerheten til parameterverdiene
v. Andre typer avgjørelser som tas på bakgrunn av modellog parameter-usikkerhet. (Risikoanalyse)
i.
Statistisk skoler- Frekventistisk
Klassisk/frekventistisk: Kun data tilordnes en
sannsynlighetsfordeling. Ofte basert på likelihood, f(D|).
Fokus på estimering ved kun å bruke data og modell.
Modellvalg og usikkerhetsanslag fra sannsynligheten for å
reprodusere noe som ligner på de data man fikk.
Mens parameterne selv ikke kan ha sannsynlighetsfordeling,
kan man tilordne en til estimatorer. En estimator er en
metode for å lage et parameter-estimat fra data. Før
dataene kommer, vil dermed en estimator ha en
sannsynlighetsfordeling.
Moment-estimering
Som nevnt, kan en stokastisk variabels forventingsverdi estimeres via
gjennomsnittet og variansen kan estimeres via observert
kvadratavvik. Forventing og varians er eksempler på momenter,
forventinger til polynomiske uttrykk av den stokastiske variabelen.
Hvis den stokastiske funksjonen har en spesifikk parametrisk form, kan
parametrene bestemmes ut ifra dette.
Eks:
For normalfordelingen, sett =gjennomsnitt og =roten av
kvadratavviket.
For lognormalfordelingen, sett =gjennomsnittet av logaritmiske verdier
og =roten av kvadratavviket, tilsvarende.
For GEV-fordelingen, som har tre parametre, trenges et moment til
(skjevhet).
OBS: For enkelte fordelinger kan det være ganske kompliserte
sammenhenger mellom parametre og momenter
Max likelihood (ML)-estimering
Poenget med ML-estimering er å finne de
parameterverdiene som gjør data mest mulig
sannsynlige.
Likelihood er sannsynlighetstettheten for data gitt
parameterverdiene, sett på som en funksjon av
parameterne. L()=f(D|).
ML-estimering er altså å finne  slik at L() får sin
maksimale verdi.
Siden log() er en monotont økende funksjon, vil optimering
over L() og l()=log(L()) gi samme resultater. Dette kan
være hensiktsmessig, siden uttrykkene kan være enklere
etter en log-transformasjon.
ML-optimering for
normalfordelingen
Skal her gjøre en ML-optimering for normalfordelingen
som et eksempel på dette.
Anta vi har et datasett (X1,..,Xn), slik at Xi~N(,)
uavhengig. Skal estimere  og .
n
l (  ,  )  log( f ( D |  ,  ))  log(
i 1
n
1
 log(2 )  n log( )  2
2
2
 ( xi   ) 2 
1
 
exp 
2
2

2 


n
 (x  )
2
i
For at et sett estimat av  og  skal optimere l(,),
må begge deriverte være null:
l
1
 2
 
n
 ( xi   ) 
i 1
i 1
nx  n
2
l
n 1 n
   3  ( x   )2  0

  i 1

n
2
(x  )  0
 ˆ  x
1 n
 ˆ 
( x  ˆ ) 2

n i 1
Ganske så likt med
moment-estimering!
ML-optimering når ting blir
vanskelige
Ikke alle modeller gir en likelihood som lar seg
analytisk optimere.
Da blir man avhengig av å kjøre en numerisk
optimering. Her finnes det mye rart, men det
meste kan deles i to kategorier:
1. Hill-climbing/lokal klatring: Disse metodene
starter i et punkt i parameter-rommet og
bruker den lokale ”topografien” til likelihoodfunksjonen til å finne den nærmeste toppen.
Eksempel: Newton’s algoritme, Nelder-Mead.
2.
Globale metoder: Disse er mye mer
sofistikerte/kompliserte. De trenger lang
kjøringstid og ofte mye finjustering.
Eksempel: simulated annealing, genetiske
algoritmer.
Parameter-usikkerhet
Et estimat er ikke sannheten. Det kan være mange mulige
parameter-verdier som er tilnærmet like rimelige, gitt de
dataene du har.
Frekventistisk statistikk opererer med konfidens-intervaller. Et
95% konfidensinterval er en lagd fra en metode for å lage
intervaller som, før data, har 95% sannsynlighet for å
omslutte den riktige parameterverdien.
(Et Bayesiansk troverdighetsintervall har 95% sannsynlighet for å omslutte
riktig parameterverdi, gitt data).
Konfidensintervaller dannes gjerne ved å se på fordelingen til
estimatorene.
Parameter-usikkerhet - teknikker
•
•
•
Eksakte teknikker. Dette får man til når man eksakt kan regne ut
fordelingen til estimatorene. Eks. 95% konfidensintervall for
normalfordelingen fås som
( x  tn1 (0.975)s / n , x  tn1 (0.975)s / n )
der s er roten av kvadratavviket og tn-1 er den såkalte t-fordelingen med
n-1 frihetsgrader.
Asymptotisk teori. Når antall data
går mot uendelig, vil dette gjelde:
2
 l ( )
ˆ ~ N( , I( ) -1 ) der I ( )   E
er Fisher's informasjonsmatrise.
2

Dermed vil (ˆ 1.96sd (ˆ),ˆ  1.96sd (ˆ)) være et 95% konfidensintervall.
Bootstrap. Her forsøker man å gjenskape fordelingen man trekte fra,
enten ved å trekke data påny med tilbaketrekning eller ved å bruke
parametriske anslag og trekke fra modellen. Man ser på spredningen
av nye parameter-estimat.
Statistisk skoler- Bayesiansk
Bayesiansk statistikk: Her angir man en førkunnskap om
parameterverdiene, , og evt. også modellene, M. Dette
oppdateres så med data, D, via Bayes formel:
f ( D |  , M ) f ( | M )
for parameter- inferensgitt modell
f (D | M )
f (D | M) Pr(M)
Pr(M| D) 
for modell- inferens
f (D)
f ( | D, M ) 
Førkunnskap
Likelihood
Fra en førkunnskap + data får man en såkalt a’posteriorifordeling for parameterne gitt modell. Dette oppsummerer
all kunnskap man har om parameterne.
All inferens gjøres med sannsynlighetsberegninger. Når analytiske metoder
ikke strekker til, kjører man simuleringemetoder. Det er ofte mulig å
trekke fra a’ posteriori-fordelingen selv om det er biter av den (f(D)) som
ikke er analytisk tilgjengelig.
Bayesiansk vs frekventistisk
Bayesiansk
statistikk
Fordeler
Ulemper
Faglig kunnskap kan tas i bruk.
Siden du må oppgi en førkunnskap, tvinges
du til å lage meningsfulle modeller.
Resultatene er ofte lett å forstå og henger
sammen med dagligdags bruk av
sannsynlighet.
Svært kompliserte modeller kan bygges og
analyseres.
Du trenger ikke ta stilling til om noe er
fundamentalt stokastisk eller ikke.
Du får parameterusikkerhet ”gratis”.
Du blir tvunget til å oppgi en førkunnskap.
Ingen førkunnskap nødvendig, betyr en
mer ”objektiv” metode.
Frekventistisk
statistikk
Mange ferdigmetoder klare til å tas ibruk.
Med andre ord en stor ”verktøykasse”
som kan anvendes med en gang.
Enklere beregninger betyr at det er enklere
å komme i gang med bruken.
Siden førkunnskapen gjerne har en subjektiv
karakter, blir resultatet å anse som
subjektivt også.
Ofte ikke så mange ferdigmetoder
tilgjengelig.
Regningen før du får resultater er oftere
vanskelig.
Vanskelig å benytte relevant faglig
førkunnskap.
Vanskelig å forstå hva resultatene faktisk betyr!
Kompliserte modeller kan være nærmest
umulig å analysere med frekventistiske
metoder.
Du må ta stilling til om noe er fundamentalt
stokastisk eller ikke.
Parameterusikkerhet er en separat oppgave du
må gjøre etter estimering.
Frekventistisk estimering kan inneholde ”bugs”,
sett i vannføringskurve-estimering.
Bayesiansk vs frekventistisk –
det pragmatiske aspektet
Når modellkompleksiteten er under en hvis terskel, er frekventistisk
metodikk enklest. Over terskelen blir det enklere med Bayesiansk
metodikk.
Arbeid
Frekventistisk
Bayesiansk
Kompleksitet
Modell-testing
•
•
Iblandt er vi ikke sikre på hvilken modell vi skal bruke.
Hvis det er snakk om vi trenger en spesifikk parameter eller ikke (kan vi
anta at forventingen i en normalfordeling antar en spesifikk verdi, f.eks.),
kjører vi modelltesting. Dette gjøres ved:
1.
2.
Sjekk om et 95% konfidensintervall omslutter den spesifikke verdien.
Beregn en p-verdi: sannsynligheten for å få en like ekstrem verdi som den vi
fikk, gitt at nullhypotesen (den spesifikke parameterverdien) stemmer. Går
denne under 5%, forkaster vi nullhypotesen med 95% konfidens. P-verdier kan
bruke trinnløst til å angi hvor sterk/svak nullhypotesen er.
En metode for å kjøre hypotesetesting er likelihood-ratio-metoden:
2(l (ˆ)  l0 (ˆ0 )) ~  p
der p er forskjellen i antallparametre,ˆ er ML - estimatetfor
full modellog ˆ0 er ML - estimatettil redusert modell,der en
del parametreblir tilegnet spesifikkeverdier.
•
Et alternativ når vi ikke har såkalte ”nøstede” modeller er
informasjonskriterer (AIC,BIC), der  2l (ˆ)  staffeledd( p) der p  antallparametre
blir minimert.
Når modell krasjer med virkeligheten
Ønsker å lage konfidensintervall for gjennomsnittelig
mammut-masse
Datasett: x=(5000kg,6000kg,11000kg)
Modell 1: xi~N(,2) i.i.d.
 Tillater mammuter å ha negativ masse!
 Resulterer i 95% konfidens-intervall, C()=(-650kg,15300kg) inneholder
verdier som bare ikke kan stemme.
Modell 2: log(xi) ~ N(,2) u.i.f. (xi ~ logN(,2) )
 Kun positive målinger for forventninger tillatt, E(xi)=exp(+2/2).
 Resulterer i 95% bootstrappet konfidens-intervall: (2500kg,10400kg).


Enda bedre hvis vi kan legge til førkunnskap.
Å få et forventningsrett estimat er dog vanskeligere . Hvis kun dette
er ønsket, kan modell 1 være bedre.
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
sammeheng mellom repsons 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  1x1,i  2 x2,i  ...  p x p,i   i der  i ~ N (, )
er uavhengig støy.
Lineær regresjon - eksempel
vekt
 Eksempel:
vekti  a  b * høydei   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
vekti  A* høydebi * Ei der Ei ~ log N (0, )
høyde
Lineære regresjon – når man går
amok i forklarinsvariable
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å
høyere-ordens polynomledd i
høyde-mot-vekt-eksempelet:
vekt
vi  0  1hi  2hi  2hi  4hi   i
2
3
4
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 parameterusikkerhetene 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. R rapporterer
p-verdier for all forklaringsvariable.
Det siste kan gjøres ved å:
a)
b)
c)
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 å velge den med best
informasjonskriterie. (Ikke anbefalt!)
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  1x1  2 x2  ...  p x p og Yˆ  ˆ0  ˆ1x1  ˆ2 x2  ... ˆ p x p
samt usikkerheten til prediksjoner:
Yˆpred  ˆ0  ˆ1x1  ˆ2 x2  ... ˆ 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 å 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
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
En nærmere titt gir
optimal h0=+8cm.
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 2
– 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. Detter er sammensatt av
kobinasjoner av modeller som harfø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 uavhengigstøy
MA-modeller: Modeller basert på glidende midling:
X t   t  1 t 1  ...  p t  p der  t ~ N (0, ) er uavhengigstø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.
En AR-modeller kan undersøkes ved et tilsvarende plott kalt ”partial
autocorrelation function” (pacf). Data produsert av en AR-modell vil
ha et pacf plott med bare et lite antall 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.