Kurs i praktisk bruk av Bayesianske metoder.

1. Bayesiansk statistikk og simulering Kurs i praktisk bruk av Bayesianske metoder.

2. Hvorfor statistikk? • Problemstillinger med ukjente størrelser: modellparametre, tidsserie-størrelser, kalibrering av hydrologiske modeller. • Trenger metoder for å håndtere det ukjente ut ifra det kjente: målinger. • Målingene har gjerne et innslag av tilfeldigheter. Kan ikke på forhånd spå hva utfallene skal bli -> sannsynlighet. • Statistikk: Fag der man fra målinger forsøker å si noe om modellene.

3. Bayesiansk statistikk • Bayesiansk statistikk anvender sannsynligheter hele veien. Ukjente størrelser håndteres via det som er kjent og tilknyttet de ukjente størrelsene. Sannsynlighet for modell regnes ut via sanns. for måledata via Bayes formel: Pr(M | D) = Pr(D | M) Pr(M) / Pr(D) • Hvis det er modell-parametre, , det er snakk om, fås tilsvarende: f(|D) = f(D| ) f() / f(D) • Krever a) en modell som angir f(D| ). b) en førkunnskap om modellen, f(). • Får fordelingen til parameterene etter datahåndtering. Altså vår kunnskap om modellen etter at data er blitt lagt til.

4. Hvorfor Bayesiansk statistikk? • I hydrologi og geofysiske fag generelt finnes for oftest førkunnskap om modeller og modell-parametre. • Resultatene fortolkes ofte enkelt av ikke-statistikere. • Mer egnet til risikoanalyse. • Kan håndtere svært komplekse modeller

5. Bayesianske resultater • Primært resultat: f( |D) (a’ posteriori-fordelingen). Alt annet avledes av dette. • Kan lage histogram og spredningsplott av f( |D) • Kan finne fordelingen til avledede størrelser, =g(). • Kan finne fordelingen til nye målinger, prediksjon. • Kan lage et ‘beste estimat’ for  via fordelingens forventningsverdi eller median. • Kan angi spredningen via standardavvik eller troverdighetsintervaller. • Kan benytte resultatet til å beregne risiko=forventet tap.

6. Problemer med Bayesiansk statistikk • Det kan være lenger vei til målet i Bayesiansk statistikk. • Spesifisering av førkunnskap i form av en fordeling kan være vanskelig. • Normaliseringskonstanten i Bayes formel: f(D) =  f(D| )f()d, kan være umulig å regne ut analytisk. • Median, forventning, standardavvik og troverdighetsintervall til f(|D) er ikke alltid analytisk tilgjengelig.

7. Probabilistisk tenkning • Håndterer utsagn med ukjent sannhetsgehalt via sannsynligheter, f. eks. sannsynligheten for ”det blir regn i morgen” eller ”det er liv på Mars”. • Pr(utsagn)=1 angir at vi er sikker på at utsagnet er sant. Tilsvarende: Pr(utsagn)=0 angir at det er usant. Alt imellom angir ulik grad av usikkerhet. • Når nye data kommer inn, oppdateres sannsynlighetene ved å se på sannsynligheten for utsagnet betinget på det du nå vet, Pr(utsagn | data). Håndteres via basislovene for sannsynlighet.

8. Basislover for sannsynlighet 1) 0Pr(A)  1 for alle utsagn A. Sannsynlighet måles fra 0-100%. 2) Pr(A) + Pr(ikke A)=1. Sannsynlighet for at det regner pluss sannsynlighet for at det ikke regner er lik 1 (sikkert). 3) Pr(A og B) = Pr(A|B)Pr(B) = Pr(B|A)Pr(A) Lov for betinget sannsynlighet. Gir også definisjonen av uavhengige hendelser. A og B er uavhengig hvis Pr(A|B)=Pr(A) => Pr(A og B)=Pr(A)Pr(B)

9. Nyttige avledede regler 4) Pr(B1 eller B2 eller B3)= Pr(B1)+Pr(B2)+Pr(B3) hvis B1, B2 og B3 er gjensidig utelukkende. Eks: Sannsynlighet for å få 2, 4 eller 6 på terningen er Pr(toer)+Pr(firer)+Pr(sekser)=3/6=1/2. 5) Pr(A eller B)=Pr(A)+Pr(B)-Pr(A og B) generelt. Eks: Pr(partall eller over to på terningen)= Pr(partall)+Pr(over to)-Pr(partall over to)= ½+4/6-2/6=5/6. 6) Hvis B1,…,Bk er gjensidig utelukkende, og dekker alle muligheter (partisjon), så er Pr(A)=Pr(A|B1)Pr(B1)+…+Pr(A|Bk)Pr(Bk) Eks: Sannsynlighet for at du blir våt på en tur er sannsynlighet for at du blir våt gitt at det regner ganger sannsynligheten for at det regner + sannsynligheten for at du blir våt gitt at det ikke regner ganger sannsynligheten for at det ikke regner.

10. Bayes formel • Starter med basisregel 3, og antar at M er en modell og D er data: Pr(M og D)=Pr(M|D)Pr(D)=Pr(D|M)Pr(M) • Snur litt på denne, og får Bayes formel: Pr(M|D)=Pr(D|M)Pr(M)/Pr(D). • Kan bruke avledet regel 6 til å finne Pr(D). Pr(D)=Pr(D|M1)Pr(M1)+….+Pr(D|Mk)Pr(Mk) • Bayes formel nå:

11. Probabilistisk tenkning - eksempel • Anta at hvis det regner, så er det overskyet og at det ikke alltid regner og ikke alltid er overskyet eller ikke. • Du vet ikke om det regner eller ikke, men du ser ut vinduet og ser at det er overskyet. Er det nå større grunn til å tro at det regner? • Pr(regn|overskyet) = Pr(overskyet|regn)Pr(regn)/Pr(overskyet)= Pr(regn)/Pr(overskyet)>Pr(regn) • Altså, sannsynligheten for at det regner har økt med observasjonen ”det er overskyet”. • Tall-eksempel: Pr(regn)=20%. Pr(overskyet)=40%. Pr(regn|overskyet)=20%/40%=50%. Altså øker sannsynligheten fra 20% til 50%. I logikken er det å argumentere mot implikasjonspilen feil. Men ikke her.

12. Eksempel 2 – flom over gitt terskel. • Lurer på gjentaksintervallet til en flomhendelse der vannet flommer over noen gitte jorder. • Feltet er umålt, men en bonde forteller at de 50 år han har vært der, har det flommet slik kun to år. • Hva blir sannsynligheten for en slik flom per år og hva blir estimert gjentaksintervall?

13. Kontinuerlige modell-parametre • Når en modell består av masse modeller av samme type kun karakterisert ved ulike parameter-verdier, håndterer vi disse som vi håndterer modell-sannsynligheter. • Forskjellen er at de har kontinuerlige mulige verdier. • Benytter sannsynlighetstettheter heller enn sannsynligheter. • Bruker integraler hellers enn summer, f.eks. i regel 6:

14. Eksempel 2 - forts • Modell: Pr(k flomhendelser over n år | p) = • Førkunnskap: mer eller mindre ingen førkunnskap, f(p)=1 for 0p1. • p|D ~ beta(k+1,n-k+1) Beta-fordelingen

15. Eks. 2, resultat (plott) zoom

16. Eks 2 – estimering • Estimat på den ene parameteren vi har, p: • E p|D = (k+1)/(n+2) (=3/52=0.05769). . • Mode p|D = k/n (=2/50=0.04) . • Median p|D har inget analytisk uttrykk. Kan trekke 100000 ganger fra f(p|D)=beta(3,49) og beregne sample-medianen. Fikk i så tilfelle 0.05193 2.6/50.

17. Eks 2- usikkerhets-estimat • Standardavvik/varians. Analytisk resultat: var(p|D)=(k+1)(n-k+1)/(n+2)2(n+3). I vårt tilfelle: sd(p|D)=0.032 • Troverdighetsintervall, et intervall som med en gitt sannsynlighet omslutter riktig verdi. Typisk: 95% troverdiget. 1) Enten finne et nivå som gjør at det område der f(p|D)>nivå har sannsynlighetsmasse 0.95 2) Eller finne pmin,pmax slik at Pr(p<pmin|D)=0.025 og Pr(p>pmax|D)=0.025.

18. Eks 2 - troverdighetsintervall • Metode 2. Kan benytte seg av tabeller til å finne at Pr(p<0.0122|D)=0.025 og Pr(p>0.1346|D)=0.025.

19. Eks 2 - prediksjon Ser på hva fremtidige data kan bli, k2 av n2 Betinget på data heller enn den ukjente parameteren p, er ikke fremtidige data binomisk fordelt lenger, men med en spesielt fordeling som ligner på hypergeometrisk fordeling.

20. Eks 2 - prediktiv fordeling Ser på prediksjon av 100 nye år basert på data og basert på frekventistisk estimat, p=0.04.

21. Eks 2 - gjentaksintervall • Gjentaksintervall er definert som forventet antall år før gjentak av en hendelse, T=1/p. • Gitt fordelingen til p kan vi finne fordelingen til T, men egenskapene til denne kan være ganske ukjente. Benytter derfor simulasjon. Trekker igjen 100000 p|D~beta(2+1,48+1) og får derfor også 100000 sampler av T|D. • Kan så studere histogram, forventning, standardavvik og troverdighetsintervall basert på sample-estimat.

22. Eks 2 – gjentaksintervall-samples • Estimere forventing fra sample-gjennomsnitt: E T  25.5. • Sample standardavvik, sd(T)  25 • Troverdighetsintervall fra estimerte kvantiler (2.5% og 97.5 -kvantil): (7.4 – 81.7)

23. Eks 3 – Normalfordelte data • Skal finne forventet vannføring 24/6 for stasjon Nor (2.2.0). • Har et tidsserie-sett fra 1937 til 1997 med vannføringer for 24/6, n=61. • Antar at qi=log(Qi)~N(,2) u.i.f. Antar i utgangspunktet kjent standardavvik, =0.523. • Førkunnskap: Ønsker en førkunnskap for log vannføring på formen N(0, 2), altså spesifisert ved forventet nivå og usikkerhet. • Antar at spesifikt avløp ligger mellom 1 og 1000 l/s/km2 med 95% sannsynlighet. Nor har nedbørsfelt 19040km2, som gir et 95% troverdighetsintervall før data (prior), på ca. 19m3/s til 19000m3/s. For log vannføring blir dette (2.947,9.854). • For normalfordelingen er et 95% troverdighetsintervall (0 -1.96, 0+1.96), som gir 0=6.40 ,=1.76 for vårt tilfelle.

24. 0, 2 Grafisk modell: Hyper-parametre  Parametre • Naturlig konjugert: Førkunnskap har samme form som likelihooden’s funksjonsavhengning for parameterene. Gir ofte en kunnskap etter data på samme formen. • For binomisk fordelte data var beta-fordelingen den konjugerte. • For normalfordeling er normalfordelingen den konjugerte • Gir at den nye fordelingen igjen kan karakteriseres med hyper-parametre. q1,…,qn Data

25. Eks 3 - utregning • Førkunnskap: • Data: • Etter data:

26. Eks 3 – Før og etter • Har q=6.01, n=61 gir ~N(6.02, 0.0672) • Altså 0*=6.02, *=0.067. • Merk at forventningen er nesten lik q og varians nesten lik 2/n. Førkunnskapen spilte relativt liten rolle i dette tilfellet.

27. Eks 3 - forventet vannføring • Prediktiv fordeling: qnew|D~N(0*, 2+*2) • Har at Qnew=exp(qnew), som gir EQ=exp(0 *+(2+*2)/2) =472.6m3/s • (Gjennomsnittelig Q: 472.9m3/s) f(exp(0 *)|D) f(Qnew|D)

28. Eks 4 – normalfordelt med ukjent varians, 2 • Antar igjen qi~N(, 2 ), med ~N(0,2) men nå med ukjent 2. • Trenger altså en prior for 2 også. • Vanlig: 2 ~ IG(,) (gir enklere regning) • f(2)= /() (2)--1 exp(-/2)

29. Eks 4 – ukjent 2 – grafisk modell 0, 2 ,  Hyper-parametre 2  Parametre q1,…,qn Data I vårt tilfelle fra eks 3, tilpasser ,  slik at Pr(log(1.25)<<log(5)): =1.38, =0.215

30. Førkunnskap: Data: Etter data: Eks 4 - utregning som før Ingen analytisk normalisering!

31. Simulering • Simulering benyttes når analytiske metoder kommer til kort. • Bruker det til: • Beregne momenter og andre integral som ikke er direkte analytisk tilgjengelig. (Monte Carlo). • Trekke fra en fordeling med ukjent normalisering (MCMC).

32. Monte Carlo-metoder • Betegnelsen ”Monte Carlo metoder” benyttes på alt som bruker tilfeldige trekninger fra gitte fordelinger til å beregne noe. • Mest vanlig å benytte til å beregne integraler: • Eks: Vår beregning av ET=E1/p i eks 2. • Beregning av  ved å trekke fra uniform fordeling over [-1,1]x[-1,1] og se andelen som er innefor enhetssirkelen. /4andel.

33. Importance sampling • Metode som kan effektivisere Monte Carlo-metoder. • Skal fortsatt beregne Efg(x), men nå ved en alternativ forslagsfordeling, q. • Monte Carlo-estimat:

34. Fordeler med importance sampling • Økt effektivitet hvis q(x)f(x)g(x) • Er ikke avhengig av normaliseringen til f(x), siden vi har w(x)=f(x)/q(x) både over og under brøkstreken. Kan beregne momenter til f(|D)=f(D|)f()/f(D) uten å kjenne f(D). • Ulempe: Kan ha stabiliseringsproblemer i enkelte sammenhenger.

35. Trekning fra a’ posteriori-fordelingen • Hvis normaliseringkonstanten f(D) = f(D|)f()d ikke er kjent, kan vi likevel trekke fra f(|D) = f(D|)f()/f(D). • Metodene for å gjøre dette heter Markov Chain Monte Carlo-metoder (MCMC).

36. Markov-kjeder • En Markov-kjede er et sett tilfeldige variable der en trekning kun avhenger av forrige trekning: f(xi | xi-1,xi-2,…,x1)=f(xi|xi-1). • Eks: xi=xi-1+ei der ei~N(0,2) uif. (Random Walk, RW) (Muligens også en del av våre tidsserier) • Under visse omstendigheter vil kjeden konvergere mot en fast fordeling, f(xi)h(x) når i. • MCMC: Lager en Markov-kjede over parameter-rommet, , slik at f(i)f(|D).

37. Første MCMC-Metropolis (en parameter) • Lagd i 1952 i forbindelse med H-bombe-utviklingen. • Algoritme for å trekke fra f(|D), gitt unormalisert versjon av fordelingen, h() (typisk f(D|)f()): • Start med en trekning av (0)~g() der g er en forslagsfordeling (gjerne a’ priori-fordelingen). • Loop: i=1…N • Trekk new~p(new| (i-1)) der p er en symmetrisk forslagsfordeling, p(x|y)=p(y|x). Eks: normalfordeling eller uniform fordeling sentrert rundt (i-1). (RW-Metropolis) • Trekk u~U(0,1) • Hvis u<h(new)/h((i-1)) så sett (i)= new. Hvis ikke, sett (i)= (i-1). • Behold (N) som en trekning fra f(|D).

38. MCMC med mer enn en parameter • Kan lage forslagsfordelinger for enkelt-parametre, og så gjennomløpe rekka av parametre. • Kan lage forslagsfordelinger, p(new| old), for mer enn en parameter av gangen. (Gjort rett kan dette være svært effektivt, men mer jobb for statistikeren/programmereren!)

39. draw.theta=function(N) { # 1. Første trekning: mu=rnorm(1,mean(q),0.1) sigma2=1/rgamma(1,alpha,beta) # 2. Løkke: for(i in 1:N) { mu.new=rnorm(1,mu,1) # 3 u=runif(1) # 4 if(u<h(mu.new,sigma2)/h(mu,sigma2)) # 5 mu=mu.new sigma2.new=rnorm(1,sigma2,0.1) u=runif(1) if(u<h(mu,sigma2.new)/h(mu,sigma2)) sigma2=sigma2.new } c(mu,sigma2) # 6 - returner theta_N } theta=array(NA,c(2,100)) for(i in 1:100) theta[,i]=draw.theta(1000) mu=theta[1,] sigma=sqrt(theta[2,]) Eks 4 – R-kode # Antar vi allerede har lest inn data, q # Førkunnskap: alpha=1.38 beta=0.215 mu0=6.40 tau=1.76 ig=function(x, a, b) a^b/gamma(a)*x^(-a-1)*exp(-b/x) h=function(mu,sigma2) { if(sigma2>0) prod(dnorm(q, mu, sqrt(sigma2)))* dnorm(mu,mu0,tau)*ig(sigma2,alpha,beta) else 0 }

40. Resultat: • Har nå 100 trekninger fra f(=(,2)|D). • Kan estimere forventing og varians via gjennomsnitt og sample-varians. • E1/100 i (i)=6.01 • E0.59

41. Triks for økt robusthet – logaritmisk sanns. draw.theta=function(N) { # 1. Første trekning: mu=rnorm(1,mu0,tau) sigma2=1/rgamma(1,alpha,beta) # 2. Løkke: for(i in 1:N) { mu.new=rnorm(1,mu,1) # 3 u=runif(1) # 4 if(log(u)<log.h(mu.new,sigma2)-log.h(mu,sigma2)) # 5 mu=mu.new sigma2.new=rnorm(1,sigma2,0.1) u=runif(1) if(log(u)<log.h(mu,sigma2.new)-log.h(mu,sigma2)) sigma2=sigma2.new } c(mu,sigma2) # 6 - returner theta_N } theta=array(NA,c(2,100)) for(i in 1:100) theta[,i]=draw.theta(1000) mu=theta[1,] sigma=sqrt(theta[2,]) # Antar vi allerede har lest inn data, q # Førkunnskap: alpha=1.38 beta=0.215 mu0=6.40 tau=1.76 ig=function(x, a, b) a^b/gamma(a)*x^(-a-1)*exp(-b/x) log.h=function(mu,sigma2) { if(sigma2>0) sum(log(dnorm(q, mu, sqrt(sigma2))))+ log(dnorm(mu,mu0,tau))+ log(ig(sigma2,alpha,beta)) else -1e+200 }

42. Effektivitet – burnin og effektiv uavhengighet • Har at en MCMC konvergerer mot riktig fordeling, men hvor mange trekninger trenger vi? (burnin) • Når vi har nådde en stabil fordeling, burde vi slippe å starte på ny for å trekke en gang til. Hvis vi kan anta at (i+k) ca. uavhengig av (i), bør det holde å fortsette kjeden og beholde hver k’te trekning.

43. Burnin • Kan undersøke burnin ved å se på grafen {i,(i)}. Typisk vil kjeden i starten bevege seg mot senteret i fordelingen, før den stabiliserer seg. • Mer avansert: Kan starte flere kjeder og studere når variansen innad i kjedene blir større enn variansen mellom kjedene. Kan også bruke andre mål for å sjekke om det er forskjell på kjedene. Viktig: Start fra ulike startsteder!

44. Effektiv uavhengighet • For å trekke slippe å trekke forferdelig mange ganger, ønsker vi at antall iterasjoner, k, før effektiv uavhengighet skal bli så liten som mulig. Formel basert på autoregresiv antagelse: k=/(1-) der  er estimert ett-stegs autokorrelasjon. • For liten akseptanserate og mange trekninger må foretas før en vi kan anta uavhengighet. • For lite bevegelse kan også komme av at forslagsfordelingen gir forslag som er svært nærme forrige verdi. Akseptanseraten blir høy, men kjeden beveger seg lite for hver iterasjon. • Må derfor finne et kompromiss mellom for høy og for lav akseptanserate. Anbefalt: akseptanserate mellom 0.25 og 0.5.

45. Typisk rutine: j=j+1 k=k+1 if(j>burnin && k>=spacing) { theta[,i]=c(mu,sigma2) i=i+1 k=0 } } theta # 6 - returner theta } draw.theta=function(N, rw.mu, rw.sigma2, burnin, spacing) { # 1. Første trekning: mu=rnorm(1,mu0,tau) sigma2=1/rgamma(1,alpha,beta) theta=array(NA,c(2,N)) j=i=k=1 # 2. Løkke: while(i<=N) { mu.new=rnorm(1,mu,rw.mu) # 3 u=runif(1) # 4 if(log(u)<log.h(mu.new,sigma2)-log.h(mu,sigma2)) # 5 mu=mu.new sigma2.new=rnorm(1,sigma2,rw.sigma2) u=runif(1) if(log(u)<log.h(mu,sigma2.new)-log.h(mu,sigma2)) sigma2=sigma2.new

46. Metropolis-Hastings • I Metropolis-algoritmen er det et krav om symmetriske forslagsfordelinger. Kan fjerne dette kravet i Metropolis-Hastings. Antar vi igjen ønsker å trekke fra den unormaliserte fordelingen h(). • Start med en trekning av (0)~g() der g er en forslagsfordeling. • Loop: i=1…N • Trekk new~p(new| (i-1)) • Trekk u~U(0,1) • Hvis u<h(new)g((i-1)|new)/h((i-1))g(new|(i-1)) så sett (i)= new. Hvis ikke, sett (i)= (i-1).

47. Metropolis-Hastings eksempel: independence sampler • Anta vi vet at |D er ca. fordelte som en kjent fordeling p(). • Bruker p som forslagsfordeling, og forslår dermed ny  uavhengig av forrige trekning. • Godtar nytt forslag hvis u<h(new)g((i-1))/h((i-1))g(new) • Blir mer og mer effektiv jo høyere akseptanseraten blir. • Kan typisk finne en normalfordeling som ligner a’ posteriorifordelingen ved å trekke noen ganger via RW Metropolis og finne forventning og varians derifra.

48. Gibbs-sampling • En spesialversjon av MH-algoritmen som sørger for at akseptanseraten er lik 1. • Kun meningsfylt når man har flere parametre. • Start med en trekning av (0)~g() der g er en forslagsfordeling. • Loop: i=1…N • Loop: j=1…k der k er antall parametre • Trekk (i)j ~ f(j| D,(i)1,…, (i)j-1, (i-1)j+1,…, (i-1)k)= f(j| D,(i)-j)

49. Gibbs-sampling – eks 4 • Har at • Ser vi på funksjonsavhengigheten for  når 2er kjent, har dette formen til en normalfordeling, som vi kjenner normaliseringen til. Dette ble allerede gjort i eks 3: • Tilsvarende fås for 2når  er kjent:

50. Gibbs-sampling - kode draw.theta=function(N, rw.mu, rw.sigma2, burnin, spacing) { # 1. Første trekning: mu=rnorm(1,mu0,tau) sigma2=1/rgamma(1,alpha,beta) n=length(q) meanq=mean(q) theta=array(NA,c(2,N)) j=i=k=1 while(i<=N) { mu=rnorm(1,(mu0*sigma2/n+meanq*tau^2)/(sigma2/n+tau^2),sqrt(tau^2*sigma2/n/(sigma2/n+tau^2))) sigma2=1/rgamma(1,alpha+n/2,beta+sum((q-mu)^2)/2) j=j+1 k=k+1 if(j>burnin && k>=spacing) { theta[,i]=c(mu,sigma2) i=i+1 k=0 } } theta # 6 - returner theta }