Fasit

Report
Oppgave 1: Utfør R-koden på
http://folk.uio.no/trondr/nvekurs2/hoelen1.R
Denne koden skal gi svar på følgende:
a) Ta en titt på årsvannføringer (snitt) fra Hølen.
b) Se på histogram sammen med en normalfordeling med
samme snitt (forventing) og varians som data. Se om
dataene ser noenlunde normalfordelt ut.
Ser ok ut.
Oppgave 1: forts.
c) Se på histogram sammen med en normalfordeling med
samme snitt (forventing) og varians som data. Se om
dataene ser noenlunde normalfordelt ut.
Ser ok ut.
Oppgave 1: forts.
d) Gjør det samme som i b og d, men bruk lognormalfordelingen i stedet, der
log-snitt og log-varians er den samme fordelingen som i data.
Ser også ok ut, men kanskje ikke like bra for lave verdier (det finnes noen verdier som er ekstremt
lave ifølge lognormal-fordelingen).
e) Gjenta b-d for døgnvannføring også (finnes på
http://folk.uio.no/trondr/nvekurs2/TrendDognHoelen.txt). Hvis
konklusjonene blir litt ulike, hva er grunnen?
Ekstremt
dårlig!
Bedre, men ikke bra.
Oppgave 2: Forventingsverdien til
årsvannføringer fra Hølen.
http://folk.uio.no/trondr/nvekurs2/hoelen2.R
a) Estimer forventningsverdien.
11.67041
b) Sjekk om forventingen er 10m3/s ved en t-test
(tar hensyn til usikkerheten i estimert varians).
Bruk gjerne 5% signifikansnivå (konfidens 95%).
t = 5.4363, df = 84, p-value = 5.233e-07
alternative hypothesis: true mean is not equal to 10
95 percent confidence interval:
11.05937 12.28145
sample estimates:
mean of x
11.67041
95% konf. Int. fra 11.06 til 12.28,
omslutter ikke 10, altså kan vi
forkaste forventing=10.
p-verdi=0.5ppm, så dermed
kunne vi tatt i med mye mer
konfidens enn som så, også.
Oppgave 2 forts.
c) Vis data sammen med konfidensintervallet.
Er det en bekymring at såpass masse års-snitt ligger utenfor
konfidensintervallet?
Nei. Dette konfidensbåndet antyder vår usikkerhet angående forventingen
(langtids-snittet), ikke angående enkelt-utfall.
Er det 95% sannsynlighet for at egentlig forventingsverdi ligger innefor det
spesifikke konfidensintervallet?
Nei. Et 95% konfidensintervall kan anses som en metode som før data har
95% sannsynlighet for å omslutte egentlig verdi, hva nå enn den verdien
er. Plugger man inn data, har man ingenting å ta sannsynlighet over lenger.
Frekventistisk metodikk har ikke noe konsept om sannsynlighet for
parameter verdier. Den har sannsynligheter for estimatorer, som er
funksjoner av de data man får. Men plugger man i data i en estimator, har
man dermed sikekrhet angående estimatoren.
Oppgave 2 forts.
d) Kunne vi gjort a-c for døgndata også?
5a kan gjøres trivielt og vil uansett være et estimat på forventingsverdien.
Analysen i oppgave 5b og 5c avhenger av uavhengigshetsantagelser, som
slettes ikke kan stemme for døgndata!
e) Skal nå foreta samme analyse der vi bruker lognormal-fordelingen hellers
enn normalfordelingen. Kjør en bootstrap-analyse som angir 95%
konfidensintervall. Hva sier dette om antagelsen forventing=10m3/s?
Anslått 95% konfidensintervall omslutter ikke 10. (Det må sies at bootstrap-estimatet
av 95% konfidensintervall er litt naivt her, for å holde ting enkelt.) Konfidensintervallet
er ganske like du får med t-testen (under normalfordelingsantagelser).
Konklusjon: forventing=10m3/s stemmer nok ikke.
Ekstraoppgave 1: Terningsutfall
På en kubisk terning er det 1/6 sannsynlighet for hver type utfall fra 1 til 6. Ved to
terninger, er utfallene antatt uavhengig.
a)
Hva er sannsynligheten for å få et spesifikt utfall på to terninger, f.eks.
sannsynligheten for å få 5 på første terning og 2 på andre?
Uavhengige utfall:
P(femmer på første og toer på andre terning) =
P(femmer på første)*P(toer på andre) = 1/6*1/6 = 1/36
b) Hva blir da sannsynligheten for å få sum=2 på de to terningene? Gjenta for sum=3,
sum=4, sum=5, sum=6 , sum=7, sum=8.
Sum=2. En måte å få til det på; 1-1. Altså P(sum=2)=1/36
Sum=3. To måter å få til det på; 1-2, 2-1.
Altså P(sum=3)=P(1-2)+P(2-1)=2/36=1/18
Sum=4. Tre måter å få til det på; 1-3, 2-2, 3-1.
Altså, P(sum=4)=1/36+1/36+1/36=3/36=1/12
Sum=5. Fire måter å få til det på; 1-4, 2-3, 3-2, 4-1.
Altså, P(sum=5)=4/36=1/9
Sum=6. Fem måter; 1-5, 2-4, 3-3, 4-2, 5-1. P(sum=6)=5/36
Sum=7. Seks måter; 1-6, 2-5, 3-4, 4-3, 5-2, 6-1. P(sum=7)=6/36=1/6
Sum=8. Fem måter; 2-6, 3-5, 4-4, 5-3, 6-2. P(sum=8)=5/36
Ekstraoppgave 1 (forts): Terningsutfall
På en kubisk terning er det 1/6 sannsynlighet for
hver type utfall fra 1 til 6. Ved to terninger, er
utfallene antatt uavhengig.
c) Hva er sannsynligheten for å få sum<=4?
P(sum<=4) =
P(sum=2 eller sum=3 eller sum=4) =
P(sum=2)+P(sum=3)+P(sum=4) =
1/36+2/36+3/36=6/36=1/6
d) Hva er sannsynligheten for to like?
Seks ulike utfall med to like: 1-1, 2-2, 3-3, 4-4, 5-5, 6-6
P(to like)=P(1-1)+P(2-2)+…+P(6-6)=6/36=1/6
Ekstraoppgave 1 (forts): Terningsutfall
På en kubisk terning er det 1/6 sannsynlighet for hver
type utfall fra 1 til 6. Ved to terninger, er utfallene
antatt uavhengig.
e) Hva er sannsynligheten for å få to like og sum<=4?
P(to like og sum<=4) = P(1-1)+P(2-2)=2/36=1/18
f) Hva er sannsynligheten for enten å få sum<=4 eller to like
terninger? Du kan bruke svaret fra c, d og e.
1). P(to like eller sum<=4) =
P(to like)+P(sum<=4)-P(to like og sum<=4) =
1/6+1/6-1/18 = 5/18
2). P(to like eller sum<=4)=
P(1-1)+P(1-2)+P(2-1)+P(2-2)+P(1-3)+P(3-1)+P(3-3)+P(4-4)+
P(5-5)+P(6-6)=10/36=5/18
Ekstraoppgave 1 (forts): Terningsutfall
På en kubisk terning er det 1/6 sannsynlighet for hver type utfall
fra 1 til 6. Ved to terninger, er utfallene antatt uavhengig.
g) Både fra regelen for betinget sannsynlighet og fra listen av utfall der
sum<=4, hva blir sannsynligheten for to like gitt sum<=4?
1). P(to like | sum<=4)=
P(to like og sum<=4)/P(sum<=4)=
(1/18) / (1/6) = 1/3
2). Liste av utfall der sum<=4: 1-1, 1-2, 1-3, 2-1, 2-2, 3-1
Hvert utfall har lik sannsynlighet, 1/6. To utfall med to like, 1-1 og 2-2.
Derfor P(to like | sum<=4)=2/6=1/3
h) Regn ut sannsynligheten for sum<=4 gitt to like, både fra liste av
mulige utfall og fra Bayes formel.
1). Liste over utfall med to like: 1-1, 2-2, 3-3, 4-4, 5-5, 6-6. Seks utfall
med like stor sannsynlighet. To av disse har sum<=4. Ergo
P(sum<=4 | to like)=1/3
2). Bayes formel:
P(tolike | sum  4)P(sum 4) (1 / 3) * (1 / 6)
P(sum  4 | to like) 

 1/ 3
P(tolike)
1/ 6
Ekstraoppgave 2:
På Blindern er det slik at det er 33.9% sjanse for at det
regner en dag, hvis det regnet gårsdagen, og 12.9% sjanse
for at det regner en dag hvis det ikke regnet gårsdagen. PS:
Antar stasjonaritet, altså at alle sannsynligheter er de samme uavhengig av
tidspunkt, under de samme forutsetningene.
a) Hva er sannsynligheten for at det regner en tilfeldig dag? (I.e.
hva er marginalsannsynligheten for regn?) Tips…
P(regn i dag)=
P(regn i dag og regn i går)+P(regn i dag men ikke i går)=
P(regn i dag|regn i går)*P(regn i går) +
P(regn i dag | ikke regn i går)=
33.9%*P(regn i går)+12.9%*(1-P(regn i går)= (stasjonaritet)
33.9%*P(regn i dag)+12.9%*(1-P(regn i dag)
P(regn i dag)*(100%-33.9%+12.9%)=12.9%
P(regn i dag)=12.9%/79%=16.3%
Ekstraoppgave 2 forts.
På Blindern er det slik at det er 33.9% sjanse for at det
regner en dag, hvis det regnet gårsdagen, og 12.9%
sjanse for at det regner en dag hvis det ikke regnet
gårsdagen. PS: Antar stasjonaritet, altså at alle sannsynligheter er de
samme uavhengig av tidspunkt, under de samme forutsetningene.
b) Hvorfor er sjansen for at det regnet i går gitt at det regner i
dag også 33.9%? (Tips: Bayes formel)
P(regni går | regn i dag) 
P(regni dag | regn i går) * P(regni går)

P(regni dag)
P(regni dag | regn i går) * P(regni dag)

P(regni dag)
P(regni dag | regn i går)  33.9%
(stasjonaritet)
Ekstraoppgave 3 – betingete
sannsynligheter
Hobbitun-rådet har avgjort at man skal
ekspandere hobbit-landen vestover.
Dessverre viser det seg at landene vestover
er infisert av drager!
Av de 10kmx10km arealene som er studert så
langt, var 70% av dem drage-infisert.
En standard-protokoll for område-undersøkelse
ble lagt. Et standardisert testområde av
mindre størrelse, inne i området man
undersøker, blir finkjemmet av feltbiologer.
Hobbitun biologiske avdeling har funnet at
sannsynligheten for å finne drager i et
testområde hvis området det er i er infisert
av drager, er 50%
Hvis det ikke er noen drage i området, blir det
selvfølgelig ingen deteksjon i testområdet.
Hobbit
Dragon
No dragons
Here be dragons
?
?
Ekstraoppgave 3 forts.
Modell: Områdets drage-status (L)  Sanns. for drage detektert i testområde (D)
Hva er (marginal) sannsynlighet for å
finne en drage, hvis du ikke vet om
området er infisert eller ikke?
(Hint: Loven om total sannsynlighet)
eller
Loven om total sannsynlighet:
P(finne drage)=P(finne drage | drageinfisert område)*P(drageinfisert område) +
P(finne drage | ikke drageinfisert område)*P(ikke drageinfisert område) =
50%*70%+0%*30% = 35%
Ekstraoppgave 3 forts.
Vis med Bayes formel at
sannsynligheten for å at et
område er infisert av drager, gitt
at du fant en drage i testområdet,
er 100%.
Drager i
området
Drager
funnet
Ingen
drager
Drager
P (drageinfisert område| fant drage) 
P (fantdrage | drageinfisert område)* P (drageinfisert område)

P (fantdrage)
50%* 70%/35% 100%
Finn sannsynligheten for at det er
Drager i
drager i området gitt at du ikke
området
fant noen. Kunne du forvente at
Ingen
Drager
sannsynligheten minsket fra
drager
funnet
originalsannsynligheten (70%)
selv uten å vite deteksjonssannsynligheten?
P(drageinfisert område| fant ikke drage) 
P(fantikke drage | drageinfisert område)* P(drageinfisert område)

P(fantikke drage)
50% * 70% /(1  35%)  35% / 65%  7 / 13  53.85%
Drager i
området
Ingen
drager
Siste spm, hint: Se grafisk
eller bruk evidensreglene.
Oppgave 3: Uavhengighet, Markov-kjeder og
nedbør på Blindern
Skal sjekke om tersklet døgnnedbøren på Blindern er uavhengig eller ikke.
Alternativet er at den er en Markov-kjede. Angir nedbørstilstanden (0 eller 1)
med xi, der i angir dagen. Siden vi trenger å betinge på foregående tilstand
setter vi til side den første tilstanden som x0. Antall dager bortsett fra dette
betegnes n.
Null-hypotese: Nedbørstilstanden (0 eller 1) er uavhengig fra dag til dag. Vi kan
derfor spesifisere modellen med en parameter, p=sannsynligheten for regn en
dag. Sannsynligheten for en gitt kombinasjon slike tilstander er da
P( x1,, xn )  p k (1  p)nk der k=antall dager med regn. Dette blir altså
likelihood’en til denne modellen. ML-estimatet til p blir da (kan også tas
direkte fra store talls lov): pˆ  k / n
Alternativ hypotese: Nedbørstilstanden (0 eller 1) er en Markov-kjede. Har en
overgangssannsynlighet fra ikke nedbør til nedbør på pNR og en sannsynlighet
for å at det regner neste dag hvis det regner nå på pRR. (De to andre
sannsynligheten er bare negasjonen av de foregående). Sannsynligheten for en1-pRR
gitt kombinasjon av tilstander er da:
kN
kR
P( x1,, xn )  pRR
(1  pRR )nR kR pNR
(1  pNR )nN kN
R
R
NR
R
pNR
N
der kR og nR er antall regndager og dager totalt der det regnet foregående
dag, og kN og nN er tilsvarende der det ikke regnet foregående dag. MLestimatet på de to parametrene blir da pˆ  k / n , pˆ  k / n
RR
pRR
N
N
1-pNR
Oppgave 3: Uavhengighet, Markov-kjeder
og nedbør på Blindern (forts)
Skal sjekke om tersklet døgnnedbøren på Blindern er uavhengig eller
ikke. Alternativet er at den er en Markov-kjede. Angir
nedbørstilstanden (0 eller 1) med xi, der i angir dagen. Siden vi trenger
å betinge på foregående tilstand setter vi til side den første tilstanden
som x0. Antall dager bortsett fra dette betegnes n.
a) For Blindern-data, hva blir estimatet på p=16.3% (null-hypotese)
og pNR=12.9 og pRR=33.9% (alternativ hypotese).
pRR
b) p’=16.3% (litt forskjell videre ut i desimalene).
c) 2*(la-l0)=138. 2(0.95)=3.84 => Forkaster null-hypotese med 95%
konfidens. P-verdi=7.6*10-32. Med andre ord, nei, regntilstanden er
ikke uavhengig fra dag til dag!
R
d) 95% konfidensintervall for pRR: (30.1%,37.7%)
95% konfidensintervall for pRR: (11.7%,14.1%)
1-pRR
pNR
e) For p, teoretisk: (15.1%, 17.5%)
N
For p, bootstrappet: (15.1%, 17.5%)
For p’, bootstrappet: (14.8%, 17.8%) . Rimelig at estimatet er mer
usikkert, modellen har flere parametre og siden modellen antyder at
1-pNR
det er avhengighet i data.
Oppgave 4: Medisinsk eksempel oversatt til språkbruken i Bayesiansk statistikk. (Her må man
oversette sannsynlighetstettheter tilbake til sannsynligheter og integral til summer.) Det er 0.1%
av befolkningen som har en gitt sykdom. En test gir alltid positivt utslag hvis du har sykdommen
men kun 1% sannsynlighet for positivt utslag hvis du ikke har den.
a) Hva er a’ priori-fordelingen?
Det vi ønsker å kjøre inferens på er om man er frisk eller syk. Førkunnskapen vår er at det er 0.1%
sannsynlighet for å være syk. P(syk)=0.1%, P(frisk)=99.9%.
b) Hva er likelihood’en for ulike utfall?
Mulige utfall: positiv test (pos) eller negativ test (neg). Likelihood er sannsynligheten for utfallene
(data) gitt det vi skal kjøre inferens på, altså sykdomsbildet. P(pos|syk)=100%, P(neg|syk)=0%,
P(pos|frisk)=1%, P(neg|frisk)=99%.
c) Beskriv og beregn (via loven om total sannsynlighet) a’ priori prediksjonsfordeling
(marginalfordeling).
Marginalfordelingen er sannsynligheten for det vi kjører inferens over (sykdomsstatus) ubetinget på
data-utfall, altså det vi forventer kun gitt førkunnskapen. Loven om total sannsynlighet:
P(pos)=P(pos|syk)P(syk)+P(pos|frisk)P(frisk)=100%*0.1%+1%*99.9%=1.099%
P(neg)=P(neg|syk)P(syk)+P(neg|frisk)P(frisk)=0%*0.1%+99%*99.9%=98.901%.
Alt: P(neg)=100%-P(pos))=98.901%.
d) Hva blir a’ posteriori-fordelingen, gitt positiv test eller negativ test?
Positiv test: P(syk|pos)=P(pos|syk)P(syk)/P(pos)=100%*0.1%/1.099%=9.099%
P(frisk|pos)=100%-P(syk|pos)=90.9%
Negativ test: P(syk|neg)=P(neg|syk)P(syk)/P(neg)=0% P(frisk|neg)=100%-P(syk|neg)=100%
Oppgave 4: Medisinsk eksempel oversatt til språkbruken i Bayesiansk statistikk. (Her må
man oversette sannsynlighetstettheter tilbake til sannsynligheter og integral til
summer.) Det er 0.1% av befolkningen som har en gitt sykdom. En test gir alltid
positivt utslag hvis du har sykdommen men kun 1% sannsynlighet for positivt utslag
hvis du ikke har den.
e) Øvelse i avledet størrelse (risikoanalyse): Samfunnskostnaden (K) ved en helbredende
operasjon er 10.000kr, mens nytte-effekten (N) er 100.000kr hvis man er syk, 0kr
hvis man er frisk. Hvis man har testet positiv, er det da samfunnssnyttig å starte en
operasjon med en gang? Beregn forventet samfunnsnytte, E(N-K|positiv test).
K=10.000kr, N=100.000kr hvis syk, N=0kr hvis frisk. E(N-K|pos)=(N-K|syk)P(syk|pos)+(NK|frisk)P(frisk|pos)= N*P(syk|pos)-K=
10.000kr+100.000kr*P(syk|pos)=-10.000kr+100.000kr*9.099%= -901kr. Altså
lønner ikke operasjonen seg nå, i forventning.
f) Den første testen var positiv. Hva blir a’ posteriori prediksjonsfordeling for en ny test
nå?
P(ny pos|pos)=P(ny pos|syk,pos)*P(syk|pos)+P(ny pos|frisk,pos)*
P(frisk|pos)=100%*9.099%+1%*90.9%=10.001%
P(ny neg|pos)=1-P(ny pos|pos)=90%
g) En ny test foretas og gir også positivt utfall. Hva blir nå a’ posteriori sannsynlighet for
at man er syk? Lønner det seg nå å foreta operasjonen?
P(syk|ny pos, pos)=P(ny pos | syk,pos)P(syk|pos)/
P(ny pos|pos)=100%*9.099%/10.001%=90.9%
E(N-K|ny pos,pos)=(N|syk)*P(syk|ny pos,pos)-K=100.000kr*90.9%-10.000kr=80.917kr
En operasjon vil nå i forventning lønne seg.
Oppgave 5: Forveningsverdien til årsvannføringer fra Hølen – Bayesiansk
analyse
http://folk.uio.no/trondr/nvekurs2/hoelen3.R
Antar at data er normalfordelt. Har en vag men informativ prior for
vannførings-forventningen, 0==10, se slide 17-18. Antar vi kjenner
=2.83.
a) Hvordan blir a’ posteriorifordelingen i dette tilfelle? Estimer vannførings-forventningen fra
dette. Er dette veldig forskjellig fra det du fikk i oppgave 5a?
mu.D
[1] 11.66884
11.669 nå og 11.670 tidligere. Ikke akkurat en kjempeforskjell.
tau.D
[1] 0.3068121
b) Lag et 95% troverdighetsintervall for vannførings-forventningen (Tips: 95% av
sannsynlighetsmassen befinner seg innenfor +/-1.96 standardavvik fra forventningsverdien
i en normalfordeling). Ble dette mye forskjellig fra 5b? Kan du fra dette konkludere noe
angående antagelsen vannførings-forventning=10m3/s?
c(mu.D-1.96*tau.D,mu.D+1.96*tau.D)
[1] 11.06749 12.27019
Nå, 11.07-12.27, før 11.06-12.28. ingen stor forskjell. (Det at intervallet var litt bredere for ttesten skyldes nok at man ikke kjørte noen antagelser om støy-størrelsen.)
Kan strengt tatt ikke konkludere noe angående vannførings-forventning=10m3/s. Det er ikke noe
en-til-en-forhold mellom Bayesianske troverdighetsintervall og modell-testing.
Oppgave 5 – forts.
c) Skal nå teste antagelsen vannførings-forventning=10m3/s Bayesiansk.
Sammenlign marginalsannsynlighetstettheten for de data vi fikk vs
sannsynlighetstettheten når =10. Hva antyder dette?
> d1=dnorm(mean(Q),mu.0,sqrt(tau^2+sigma^2/n))
> d1
[1] 0.03932351
>
> # Sannsynlighetstetthet hvis vi vet mu=10
> d0=dnorm(mean(Q),10,sigma/sqrt(n))
> d0
[1] 4.822864e-07
d1 = marginalsannsynlighet under vår førantagelse
d0 = marginalsannsynlighet
d1>>d0 antyder at vi har fått evidens for vår førkunnskaps-modell hellers enn
antagelsen mu=10.
Oppgave 5 – forts.
d) Skal nå bruke resultatet fra c til å regne på modellsannsynligheter. Modell 0
har =10 mens modell 1 er slik som spesifisert ovenfor. Bruk
f ( D | M ) Pr(M )
Pr(M | D) 
 f ( D | M ' ) Pr(M ' )
og anta at a’ priori-sannsynligheten for hver modell er 50%. Hva blir
konklusjonen?
> p0.D=d0*p0/(d0*p0+d1*p1)
> p1.D=d1*p1/(d0*p0+d1*p1)
> c(p0.D,p1.D)
[1] 1.226443e-05 9.999877e-01
99.9988% sannsynlighet for modell 1 (førkunnskapsmodellen) vs modell 0
(mu=10). Veldig sannsynlig at mu10, med andre ord.
Oppgave 5 – forts.
e) Lag et plott over marginalfordelingen gitt ulike utfall og sammenlign med
sannsynlighetstettheten nå =10 (likelihood under modell 0). Hva sier dette om
hvilke utfall som ville være evidens for modell 0 og 1?
> c(min(x[d0>d1]),max(x[d0>d1]))
[1] 9.19 10.81
Hvis vi fikk utfall (gjennomsnitt) mellom 9.19 og 10.81, er
marginalsannsynligheten større for modell 0 enn for modell 1. Altså ville et slikt
utfall være evidens for modell null (mu=10), ellers for modell 1.
Oppgave 6: Skal se på faren for å overgå en spesifikk vannførings-verdi.
Stasjonen Gryta har hatt vannføring>1.5m3/s y=27 ganger i løpet av t=44 år.
Antar slike hendelser foregår uavhengig i tid. Altså at antall hendelser innefor en
tidsperiode er Poisson-fordelt. Bruker gjentaks-intervall, T, som parameter i
denne fordelingen. Når man har data over en tidsperiode t, vil sannsynligheten
for y hendelser bli:
(t / T ) y t / T
L(T )  P( y hendelseri løpetav tid t | T ) 
e
y!
Denne likelihood’en maksimeres med Tˆ  t
y
. (For de med analytisk
optimerings-erfaring, vis dette).
l (T )  y log(t )  y log(T )  log( y!)  t / T
l (t )
y t
l (t )
t
  2
 0  yTˆ  t  Tˆ 
T
T T
y
Tˆ
 2l (t )
y 2t

 3
2
2
T
T
T
Ey  t / T
 2l (t )
Ey 2t
I(T ) -E


 3
2
2
T
T
T
t
 I(T ) 3
T
Viser samtidig informasjons-”matrisen”.
Oppgave 6: Skal se på faren for å overgå en spesifikk vannførings-verdi.
Stasjonen Gryta har hatt vannføring>1.5m3/s y=27 ganger i løpet av t=44 år.
Antar slike hendelser foregår uavhengig i tid. Altså at antall hendelser innefor en
tidsperiode er Poisson-fordelt. Bruker gjentaks-intervall, T, som parameter i
denne fordelingen. Når man har data over en tidsperiode t, vil sannsynligheten
for y hendelser bli:
(t / T ) y t / T
L(T )  P( y hendelseri løpetav tid t | T ) 
e
y!
Denne likelihood’en maksimeres med Tˆ  t
y
.
a) Regn ut ML-estimatet: 44år/27=1.63år.
b) 95% konfidensintervall. Sd(T)=1/I(T)=T3/2/  t=t/y3/2=0.31år.
1.63år+/-1.96*0.31år = (1.01år,2.24år).
Oppgave 7: Bayesiansk gjentaksanalyse for bestemt nivå i
kontinuerlig tid.
Skal se på faren for å overgå en spesifikk vannførings-verdi. Antar slike hendelser
foregår uavhengig i tid. Altså at antall hendelser innefor en tidsperiode er Poissonfordelt. Bruker gjentaks-intervall, T, som parameter i denne fordelingen. Da får vi
(t / T ) y t / T
P( y hendelseri løpetav tid t | T ) 
e
y!
Antar invers-gamma-fordeling (siden det er matematisk behagelig å gjøre det) for
gjentaksintervallet
   1  / T
f (T ) 
T e
( )
Får da at marginalfordelingen blir:
  y  1 y
t



P( y hendelseri løpetav tid t )  
p (1  p) der p 

y 
t

(dette er den såkalte negativ binomiske fordelingen).
Oppgave 7 (forts.): Kode finnes på
http://folk.uio.no/trondr/nvekurs2/gryta_ekstrem.R


(t / T ) y t / T
f (T ) 
T  1e / T
P( y hendelseri løpetav tid t | T ) 
e
( )
y!
  y  1 y
t
 p (1  p) der p 
P( y hendelseri løpetav tid t )  
y 
t

Stasjonen Gryta har hatt vannføring>1.5m3/s y=27 ganger i løpet av t=44 år.
a) Plott a’ priori-fordeling og marginalfordeling hvis du bruker ==1 som førkunnskap.
Ser at a’prior’en topper seg rundt T=1, men er ganske vid. Derfor svært vid
marginalfordeling også. y=27 ikke spesielt oppsiktsvekkende
Oppgave 7 (forts.): Kode finnes på
http://folk.uio.no/trondr/nvekurs2/gryta_ekstrem.R


(t / T ) y t / T
f (T ) 
T  1e / T
P( y hendelseri løpetav tid t | T ) 
e
( )
y!
  y  1 y
t
 p (1  p) der p 
P( y hendelseri løpetav tid t )  
y 
t

Stasjonen Gryta har hatt vannføring>1.5m3/s y=27 ganger i løpet av t=44 år.
b) Hva blir det generelle uttrykket for a’ posteriori-fordelingen til T?
(t / T ) y t / T    1   / T
e
T
e
f ( y | T ) f (T )
y!
( )
f (T | y ) 




y

1
f ( y)

 y

 p (1  p)
y 

(   t )  y   y 1 (  t ) / T
T
e
~  ( *    y,  *    t )
(  y )
Oppgave 7 (forts.): Kode finnes på
http://folk.uio.no/trondr/nvekurs2/gryta_ekstrem.R


(t / T ) y t / T
f (T ) 
T  1e / T
P( y hendelseri løpetav tid t | T ) 
e
( )
y!
  y  1 y
t
 p (1  p) der p 
P( y hendelseri løpetav tid t )  
t
 3y 
Stasjonen Gryta har hatt vannføring>1.5m /s y=27 ganger i løpet av t=44 år.
b) Hva blir det generelle uttrykket for a’ posteriori-fordelingen til T? Plott den for Gryta for
==1 sammen med a’ priori-fordelingen. Forsøk også ==0.5 og til og med ==0 (ikkeinformativt) . Ble det noen stor forskjell i a’ posteriori-fordelingen? Sammenlign med
klassisk estimat: TML=t/y=1.63 år. T | y ~  ( *    y,  *    t )
Posterior
med
prior
Posterior med
ulike priorer
Innzoomet
versjon
Topper seg (modus-estimat) forholdsvis nære ML-estimatet, 1.55år. Forventing /( 1)=1.67år, median=1.63år er enda nærmere.
Oppgave 7 (forts.): Kode finnes på
http://folk.uio.no/trondr/nvekurs2/gryta_ekstrem.R


(t / T ) y t / T
f (T ) 
T  1e / T
P( y hendelseri løpetav tid t | T ) 
e
( )
y!
  y  1 y
t
 p (1  p) der p 
P( y hendelseri løpetav tid t )  
y 
t

Stasjonen Gryta har hatt vannføring>1.5m3/s y=27 ganger i løpet av t=44 år.
c) 95% troverdighetsintervall T | y ~  ( *    y,  *    t ) når ==1 brukes.
Resultat (1.15, 2.42)
Asymptotisk klassisk 95% konfidens: (1.01,2.24). Ikke veldig forskjellig, men den Bayesianske
analysen antyder at usikkerheten i T er skjevt fordelt, så høyere verdier virker rimelig.
Oppgave 7 (forts.):


(t / T ) y t / T
f
(
T
)

T  1e / T
P( y hendelseri løpetav tid t | T ) 
e
( )
y!
  y  1 y
t

 p (1  p) der p 
P( y hendelseri løpetav tid t )  
y 
t

d) Kan du finne prediksjons-fordelingen til antall nye flommer på Gryta de neste
hundre år? Plott i så tilfelle denne. Sammenlign med Poisson-fordeling hvis man tar
ML-parameteren for gitt. Hvorfor er sistnevnte fordeling skarpere enn den
Bayesianske prediksjonsfordelingen?
Siden a’ posteriorifordelingen kommer fra samme fordelingsfamilie som a’
priorifordelingen, kan vi bare bytte ut  med *=+y og  med *=  +t i
marginalfordelingen.
P( y hendelseri løpet av tid t ) 
t ny
  *  yny  1 y


 p (1  p) der p 


yny
t ny   *


Poisson med gitt parameter skarpere enn marginalfordelingen siden førstnevnte ikke tar hensyn til
parameter-usikkerheten.
Oppgave 7 (forts.)


(t / T ) y t / T
f (T ) 
T  1e / T
P( y hendelseri løpetav tid t | T ) 
e
( )
y!
  y  1 y
t



P( y hendelseri løpetav tid t )  
 p (1  p) der p  t  
y


e) Kjør en enkel MCMC-algoritme fra a’ posteriori-fordelingen. Se etter når
trekningen stabiliserer seg (burn-in) og hvor mange trekninger som trenges før du
få en trekning som er ca. uavhengig (spacing).
Burnin=20-40, spacing=10
Oppgave 7 (forts.)


(t / T ) y t / T
f (T ) 
T  1e / T
P( y hendelseri løpetav tid t | T ) 
e
( )
y!
  y  1 y
t



P( y hendelseri løpetav tid t )  
 p (1  p) der p  t  
y


f) Hent 1000 uavhengige trekninger etter burn-in. Sammenlign med teoretisk a’
posteriori-fordeling (histogram og qq-plott).
Ser bra ut på min trekning.
g) Foreta ny MCMC-trekning men bruk nå a’ priori som er
f(T)=lognormal(=0,=2). (Dette kan ikke løses analytisk). Sammenlign med de
trekningene du fikk i d.
Ganske likt, men
Forholdsvis
uavhengig
kanskje noen
forskjeller helt ute i
øvre hale.
Det er helt ok…
Oppgave 8: Ekstremverdi-analyse på Bulken (rundt 120 år med data).
Kode: http://folk.uio.no/trondr/nvekurs2/bulken_ekstrem.R
Data: : http://folk.uio.no/trondr/nvekurs2/bulken_max.txt
Skal bruke Gumbel-fordelingen som fordelings-kandidat her:
f (x | ,  ) 
a)
1

e
 ( x   ) /  e ( x ) / 
Foreta et ekstremplott, det vil si sorter vannføringene og plott dem mot
n  0.12
estimert gjentakintervall t  i  0.44 der n er antall år og i er en løpe-indeks
fra n til 1.
i
b)
Foreta en ekstremverditilpasning via første to l-momenter, 1 og 2.
Sammenlign med det du får fra DAGUT. Parameterne forholder seg til lmomentene som = 2/log(29, = 1-0.57721. Estimater for 1 og 2 fås
n
som
1 n
1
Sorterte data
ˆ1 
x
n
j 1
j
 x , ˆ2 
n(n  1)
 (( j 1)  (n  j)) x
j 1
DAGUT: God overenstemmelse med beta…
( j)
> c(mu.lmom, beta.lmom)
[1] 305.22055 67.96247
Oppgave 8 –forts.
f (x | ,  ) 
1

e
 ( x   ) /  e ( x ) / 
c) Plott flomstørrelse som funksjon av gjentaksintervall gitt l-momentestimatene sammen med data (a).
Ser ok ut.
d) Foreta ML-estimering av parameterne.
> c(mu.ml, beta.ml)
[1] 304.49483 74.13713
e) Plott flomstørrelse som funksjon av gjentaksintervall gitt ML-estimatene.
Ganske like l-moment-estimatet for
små verdier. Litt uenighet for store
gjentaksintervall.
Oppgave 8 – forts.
f (x | ,  ) 
1

e
 ( x   ) /  e ( x ) / 
f) (Obs: hvis dette blir for mystisk, slutt her.) Foreta Bayesiansk analyse med flat
prior. Foreta 1000 MCMC-trekninger (burnin=1000, spacing=1000).
Sammenlign.
Histogram angir at mu er en plass
mellom 290 og 320 og at
fordelingen er noenlunde
normalfordelt.
Beta ser ut til å være et sted
mellom 65 og 85 og ser også ut
som det kan være normalfordelt.
Det ser ut fra tidsserie-plottet som
om MCMC-kjeden har stabilisert
seg men at det er en hel del
avhengighet.
Topper seg rundt mu=305, beta=75 (modus-estimat).
Ikke ulik andre estimater. # Forventings-estimat:
> c(mean(mu.mcmc),mean(beta.mcmc))
[1] 304.81509 74.86763
# Median-estimat:
> c(median(mu.mcmc),median(beta.mcmc))
[1] 305.16144 74.93275
Oppgave 8 – forts.
f (x | ,  ) 
1

e
 ( x   ) /  e ( x ) / 
g) Bruk også prediksjonsfordelingen (altså der du tar parameterusikkerheten
med i betraktningen) til å foreta samme plott som i a, c og e.
Ganske likt ML-resultatet.
Oppgave 9: Ekstremverdi-analyse på Bulken (120 år med data).
Kode: http://folk.uio.no/trondr/nvekurs2/bulken_ekstrem2.R
Data: http://folk.uio.no/trondr/nvekurs2/bulken_max.txt
Ekstrakode: http://folk.uio.no/trondr/nvekurs2/mcmc.R
Skal bruke GEV-fordelingen som fordelings-kandidat her:
x
1
x   ( 1) /  (1 (  ))1 / 
f ( x |  ,  ,  )  (1   (
))
e


a) Foreta ML-estimering av parameterne. Start med optimeringen fra
tilfeldige start-parametre (standard-normalfordelt). Gjør dette multiple
ganger for å finne den globalt maksimale likelihood’en. Plott
ekstremverdier og tilpasning.
Estimat: mu=309, sigma=76.7, kji=-0.13
b) La t være en tids-indikator som løper fra 1:120. Tilpass via ML en modell
der hver av parametrene i tur er lineært tidsavhengig, for eksempel om
(t)=0+t.
1. mu=mu0+beta_mu*t: mu0=270, sigma=72.8, kji=-0.14, beta_mu=0.68
2. sigma=sigma0+beta_sigma*t: mu=209, sigma0=71.0, kji=-0.15,
beta_sigma=0.084
3. kji=kji_0+beta_kji*t: mu=311, sigma=75.4,kji_0=-0.093,
beta_kji=-0.0009
c) AIC:1403.497 1395.430 1405.214 1405.111. Modell 1 best.
Oppgave 9: Ekstremverdi-analyse på Bulken (120 år med data).
Kode: http://folk.uio.no/trondr/nvekurs2/bulken_ekstrem2.R
Data: http://folk.uio.no/trondr/nvekurs2/bulken_max.txt
Ekstrakode: http://folk.uio.no/trondr/nvekurs2/mcmc.R
Skal bruke GEV-fordelingen som fordelings-kandidat her:
f ( x |  , ,  ) 
1

(1   (
x

( 1) / 
))
e
(1 (
x

)) 1 / 
d) theta.bayes: 309.3206264 76.7024980 -0.1311687
theta.ml: 310.5800822 76.1272248 -0.1528986
Ganske likt!
e) Tidsavhengighet i mu. Bayesiansk analyse.
Oppgave 9: Ekstremverdi-analyse på Bulken (120 år med data).
Kode: http://folk.uio.no/trondr/nvekurs2/bulken_ekstrem2.R
Data: http://folk.uio.no/trondr/nvekurs2/bulken_max.txt
Ekstrakode: http://folk.uio.no/trondr/nvekurs2/mcmc.R
Skal bruke GEV-fordelingen som fordelings-kandidatx her:

f ( x |  , ,  ) 
1

(1   (
x

( 1) / 
))
e
(1 (

)) 1 / 
e) Tidsavhengighet i mu. Bayesiansk analyse.
summary(mu1)
Min. 1st Qu. Median Mean 3rd Qu. Max.
209.4 258.9 268.7 268.5 277.8 311.8
> summary(s1)
Min. 1st Qu. Median Mean 3rd Qu. Max.
58.75 70.45 73.80 73.90 77.12 90.16
> summary(chi1)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.24890 -0.14950 -0.11730 -0.11520 -0.08459 0.14600
> summary(beta1)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.02176 0.54360 0.68970 0.68970 0.83550 1.34100
P0=21.9% P1=78.1%
Evidens for tidsavhengighet.
f) P0=58.4% P1=41.6%. Ikke overbevisende evidens.
Fordelingen til felles parametre
ser ut til å være ganske lik.
Oppgave 10:
Skal nå kjøre ARMA-tilpasning av døgndata fra Hølen.
Kode: http://folk.uio.no/trondr/nvekurs2/hoelen_arima.R
a) Plott data
b) De-trend (fjern lineær tids-trend og sesonvariasjon).
Oppgave 10 - forts.
c) Se på autokorrelsjon (og partiell autokorrelasjon).
ACF
PACF
d) Tilpass en AR(1)-modell (PS: pacf antyder at AR(2) er bedre). Se om estimert
parameter er lik noe du så i 13c.
AR(1)-koeffisient=0.9472, ganske likt
største verdi i pacf-plottet.
Coefficients:
ar1 intercept
0.9472
0.0000
s.e. 0.0018
0.0252
Oppgave 10 - forts.
e) Lag analytiske plott av residualene. Hva sier de?
Antyder at veldig mye av tidsavhengigheten er forklart allerede, men at noe
fremdeles ikke er med i modellen.
f) Forsøk så med en ARMA(1,1)-modell. Se igjen på residualene. Hva sier de nå?
At dette begynner å se bra ut?
Oppgave 11:
Bruk av Kalman-filter til å interpolere over hull på Farstad stasjon, år 1993. Skal bruke enkleste type
Markov-kjede, tilfeldig gange eller mer nøyaktig, Wiener-prosess. Dataene er har varierende
tidsoppløsning og blir log-transformert før analysen foretas.
Oppdaterings-formelen (systemligningen) er xi | xi 1 ~ N ( xi 1 ,  ti  ti 1 )
Observasjonsligningen er yi | xi ~ N ( xi , )
a)
Plott data.
Lag kunstige hull: 1993-02-01 00:00:00 - 1993-02-03 00:00:00
1993-03-01 12:00:00 - 1993-03-15 12:00:00
1993-08-01 00:00:00 - 1993-08-01 18:00:00
1993-09-01 00:00:00 - 1993-09-02 00:00:00
1993-10-01 00:00:00 - 1993-11-01 00:00:00
Plott data i hullene sammen med resterende data
Oppgave 11:
Bruk av Kalman-filter til å interpolere over hull på Farstad stasjon, år 1993. Skal bruke enkleste type
Markov-kjede, tilfeldig gange eller mer nøyaktig, Wiener-prosess. Dataene er har varierende
tidsoppløsning og blir log-transformert før analysen foretas.
Oppdaterings-formelen (systemligningen) er xi | xi 1 ~ N ( xi 1 ,  ti  ti 1 )
Observasjonsligningen er yi | xi ~ N ( xi , )
b) Ta en titt på det implementerte Kalman-filteret og sammenlign med beskrivelsen i kurset.
Det skal være som overlapp, men merk at første-instansen må håndteres på et vis. Lagde en
førkunnskap, x1~N(0,10). Siden det kun er en dimensjon i tilstands- og observasjonsligningen er det
skalar-operasjoner, ikke vektor- og matrise-operasjoner. Merk også hva som skjer når det ikke er data på
et tidspunkt. A’ posteriori = a’ priori
c) Kjør en ML-optimering og sjekk resultatene. Hva kan sies om observasjonsparameteren ?
sigma.ml=0.0330157805 tau.ml=0.0000153492
Tau blir estimert til å være forsvinnende liten. Kunne dermed brukt som observasjonsligning at yi=xi og
håndtert dette som observasjoner direkte på en tilfeldig gange.
d) Foreta en Kalman-glatting med ML-estimerte parametre. Se på total-resultatet.
Oppgave 11:
Bruk av Kalman-filter til å interpolere over hull på Farstad stasjon, år 1993. Skal bruke enkleste type
Markov-kjede, tilfeldig gange eller mer nøyaktig, Wiener-prosess. Dataene er har varierende
tidsoppløsning og blir log-transformert før analysen foretas.
Oppdaterings-formelen (systemligningen) er xi | xi 1 ~ N ( xi 1 ,  ti  ti 1 )
Observasjonsligningen er yi | xi ~ N ( xi , )
e). Sjekk hvordan interpoleringen har fungert i hullene.
Hva slags interpolasjon er dette (på log-skala). Lineær-interpolasjon.
Hva får man ut her i tillegg til interpolasjonen? Usikkerheten til estimatet!
Oppgave 11:
Bruk av Kalman-filter til å interpolere over hull på Farstad stasjon, år 1993. Skal bruke enkleste type Markovkjede, tilfeldig gange eller mer nøyaktig, Wiener-prosess. Dataene er har varierende tidsoppløsning og blir
log-transformert før analysen foretas.
Oppdaterings-formelen (systemligningen) er xi | xi 1 ~ N ( xi 1 ,  ti  ti 1 )
Observasjonsligningen er yi | xi ~ N ( xi , )
f). Kritiser modellen og se om du skjønner hvorfor tilpasningen og usikkerhetene blir slik de blir. Om du er i det
kreative hjørnet, prøv å foreslå bedre modeller.
Med forventing lik forrige verdi, er det kanskje ikke overraskende at en interpolert verdi blir et kompromiss
mellom verdiene rundt hullet. Modellen forholder seg kun til selve vannføringene, ikke til den tidsdervierte av
de (Wiener-prosessen har ikke noen derivert). At usikkerhetene bobler ut desto mer desto lenger unna
dataene vi er, er vel ikke så urimelig. At de kan bli ubegrenset store med ubegrenset store hull er vel kanskje
mindre rimelig.
For det første er modellen slik at retningen på forandringen i fremtiden ikke avhenger av retningen tidligere. En
modell som hadde en kontinuerlig forandring i ”derivert” ville derfor være bedre og gi glattere interpolering, se
hullet i starten av februar.
Ellers vil usikkerheten øke med kvadratroten av størrelsen på hullet uansett hull-mengde. Dette fordi Wienerprosessen ikke har noen stasjonær fordeling. En modell med stasjonær fordeling (og dermed gjerne vil ha en
forventning og varians som konvergerer mot noe fast) ville nok være bedre.
Ellers kan det være sesong-avhengighet både i signal og varians ville nok være mer realistisk.
En modell som inkorporerte en enkel hydrologisk modell ville nok være enda mye bedre, men betyr at mer
sofistikert metodikk tas i bruk (extended Kalman, unscented kalman eller til og med partikkel-filter).
Oppgave 12:
Bruk av Kalman-filter til å interpolere over hull på stasjonene Etna og Hølervatn. Dataene er har ekvidistant tidsoppløsning (døgn) og blir
log-transformert før analysen foretas. Skal bruke en 2-dimensjonal AR1-prosess (så den har en stasjonær tilstand), krysskorrelert støy
(slik at kompletering blir mulig) lik autokorrelasjon og varians, men individuell forventing og sesong-variasjoner i denne
forventningen (som antas lik for de to seriene).
Oppdaterings-formelen (systemligningen) er
x i | x i 1 ~ N (a x i 1  (1  a )  (t ), )
 2
  e  S1 sin(2t / 365)  C1 cos(2t / 365) 
 og   
der  (t )  
2
  h  S1 sin(2t / 365)  C1 cos(2t / 365) 
 
2 

 2 
, )
Observasjonsligningen er yi ,e / h | xi ,e / h ~ N ( x
i ,e / h
a) Plott data. Lag kunstige hull (spesifisert i T-koden). Plott data i hullene sammen med resterende data
b)
Ta en titt på det implementerte Kalman-filteret og sammenlign med beskrivelsen i kurset.
Her er vektor- og matrise-regning nødvendig, siden tilstandsrommet er to-dimensjonalt. Observasjonsrommet *kan* var to-dimensjonalt
også, men av og til har vi bare observasjoner på ett sted. Da blir observasjonsligningen annerledes (projeserer da ut kun den ene
observasjonen som finner sted). Siden observasjonsmatrisen blir gitt en fast dimensjonalitet, må også denne få sine verdier projesert ut
når det er bare en observasjon. Foretar transformasjoner av parameterverdiene slik at input-parameterene kan anta hvilke som helst
verdier på den reelle tallaksen.
Oppgave 12:
Bruk av Kalman-filter til å interpolere over hull på stasjonene Etna og Hølervatn. Dataene er har ekvidistant tidsoppløsning (døgn) og blir
log-transformert før analysen foretas. Skal bruke en 2-dimensjonal AR1-prosess (så den har en stasjonær tilstand), krysskorrelert støy
(slik at kompletering blir mulig) lik autokorrelasjon og varians, men individuell forventing og sesong-variasjoner i denne
forventningen (som antas lik for de to seriene).
Oppdaterings-formelen (systemligningen) er
x i | x i 1 ~ N (a x i 1  (1  a )  (t ), )
 2
  e  S1 sin(2t / 365)  C1 cos(2t / 365) 
 og   
der  (t )  
2
  h  S1 sin(2t / 365)  C1 cos(2t / 365) 
 
Observasjonsligningen er
yi,e / h | xi,e / h ~ N ( xi ,e / h , )
2 

 2 
c) Kjør en ML-optimering og sjekk resultatene. Hva vil autokorrelasjonen her ha å si? Hva har krysskorrelasjonen å si?
c(a.ml,mu.e.ml,mu.h.ml,sigma.ml,rho.ml,tau.ml,S1.ml,C1.ml)
[1] 0.94775606 0.16454854 -0.03336243 0.26462117 0.96971749 0.03952556 -0.07919922 -0.15421623
Tau er slettes ikke null i dette tilfellet. Det betyr en del usikkerhet selv der det er målinger. Kan være et tegn på et forbedringspotensiale i
modellen. Sterk auto- og kryss-korrelasjon, ifølge estimatet. a^30=0.5 betyr en karakteristisk tid på omtrent en måned. Krysskorrelasjonen
betyr at komplettering er mulig og kan forventes å bli forholdsvis bra.
d) Foreta en Kalman-glatting med ML-estimerte parametre. Se på total-resultatet.
De fleste hullene er små i forhold til hele serien. Man ser likevel litt usikkerhet som
tyter ut der det er større hull.
Oppgave 12:
Bruk av Kalman-filter til å interpolere over hull på stasjonene Etna og Hølervatn. Dataene er har ekvidistant tidsoppløsning (døgn) og blir
log-transformert før analysen foretas. Skal bruke en 2-dimensjonal AR1-prosess (så den har en stasjonær tilstand), krysskorrelert støy
(slik at kompletering blir mulig) lik autokorrelasjon og varians, men individuell forventing og sesong-variasjoner i denne
forventningen (som antas lik for de to seriene).
e) Sjekk hvordan kompletteringen har fungert i hullene. Hva skjer når begge stasjonene har hull samtidig? Et par kompletteringer ser ut til
å fungere mindre bra. Hvorfor?
Etna:
Her er det hull i den Hølenvatn også.
Hølenvatn:
Kompletteringen gikk mindre bra fordi
kompletteringsserien (Etna) mangler den ekstra
toppen. Usikkerheten bommer nok fordi
modellen har noe utilstrekkelig i seg.
Her er det hull i den Etna også.
Oppgave 12:
Bruk av Kalman-filter til å interpolere over hull på stasjonene Etna og Hølervatn. Dataene er har ekvidistant tidsoppløsning (døgn) og blir
log-transformert før analysen foretas. Skal bruke en 2-dimensjonal AR1-prosess (så den har en stasjonær tilstand), krysskorrelert støy
(slik at kompletering blir mulig) lik autokorrelasjon og varians, men individuell forventing og sesong-variasjoner i denne
forventningen (som antas lik for de to seriene).
Oppdaterings-formelen (systemligningen) er
x i | x i 1 ~ N (a x i 1  (1  a )  (t ), )
 2
  e  S1 sin(2t / 365)  C1 cos(2t / 365) 
 og   
der  (t )  
2
  h  S1 sin(2t / 365)  C1 cos(2t / 365) 
 
Observasjonsligningen er
yi,e / h | xi,e / h ~ N ( xi ,e / h , )
2 

 2 
f) Sjekk om krysskorrelasjonen er null (altså at det ikke er grunn til å kompletere), via likelihood-ratio-testen.
# test med 5% signifikansnivå (95% konfidens):
c(2*(l.a-l.0), qchisq(0.95,1))
[1] 521.444714 3.841459
Likelihood-forskjellen mye større enn hva man kunne forvente under null-hypotesen. Null-hypotese forkastet med 95% konfidens (og
gjerne større også).
pchisq(2*(l.a-l.0),1,lower=F)
[,1]
[1,] 2.05223e-115
P-verdi ekstremt lav! Det er veldig veldig sikkert at det er krysskorrelasjon mellom de to seriene. Siden modellen tar tidskorrelasjon med i
betraktning, skal dette resultatet ikke være en artefakt av tidsserie-aspektet til dataene.

similar documents