kurs i praktisk bruk av bayesianske metoder l.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
Kurs i praktisk bruk av Bayesianske metoder. PowerPoint Presentation
Download Presentation
Kurs i praktisk bruk av Bayesianske metoder.

Loading in 2 Seconds...

play fullscreen
1 / 58

Kurs i praktisk bruk av Bayesianske metoder. - PowerPoint PPT Presentation


  • 245 Views
  • Uploaded on

Bayesiansk statistikk og simulering. Kurs i praktisk bruk av Bayesianske metoder. Hvorfor statistikk?. Problemstillinger med ukjente størrelser: modellparametre, tidsserie-størrelser, kalibrering av hydrologiske modeller.

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'Kurs i praktisk bruk av Bayesianske metoder.' - tiva


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
kurs i praktisk bruk av bayesianske metoder

Bayesiansk statistikk

og simulering

Kurs i praktisk bruk av Bayesianske metoder.

hvorfor statistikk
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.
bayesiansk statistikk
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.
hvorfor bayesiansk statistikk
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
bayesianske resultater
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.
problemer med bayesiansk statistikk
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.
probabilistisk tenkning
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.
basislover for sannsynlighet
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)

nyttige avledede regler
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.

bayes formel
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å:
probabilistisk tenkning eksempel
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.

eksempel 2 flom over gitt terskel
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?
kontinuerlige modell parametre
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:
eksempel 2 forts
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

eks 2 estimering
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.
eks 2 usikkerhets estimat
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.

eks 2 troverdighetsintervall
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.
eks 2 prediksjon
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.

eks 2 prediktiv fordeling
Eks 2 - prediktiv fordeling

Ser på prediksjon av 100 nye år basert på data og basert på frekventistisk estimat, p=0.04.

eks 2 gjentaksintervall
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.
eks 2 gjentaksintervall samples
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)
eks 3 normalfordelte data
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.
grafisk modell

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

eks 3 utregning
Eks 3 - utregning
  • Førkunnskap:
  • Data:
  • Etter data:
eks 3 f r og etter
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.
eks 3 forventet vannf ring
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)

eks 4 normalfordelt med ukjent varians 2
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)
eks 4 ukjent 2 grafisk modell
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

eks 4 utregning
Førkunnskap:

Data:

Etter data:

Eks 4 - utregning

som før

Ingen analytisk

normalisering!

simulering
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).
monte carlo metoder
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.

importance sampling
Importance sampling
  • Metode som kan effektivisere Monte Carlo-metoder.
  • Skal fortsatt beregne Efg(x), men nå ved en alternativ forslagsfordeling, q.
  • Monte Carlo-estimat:
fordeler med importance sampling
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.
trekning fra a posteriori fordelingen
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).
markov kjeder
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).
f rste mcmc metropolis en parameter
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).
mcmc med mer enn en parameter
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!)
eks 4 r kode

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

}

resultat
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
triks for kt robusthet logaritmisk sanns
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

}

effektivitet burnin og effektiv uavhengighet
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.
burnin
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!
effektiv uavhengighet
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.
typisk rutine
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

metropolis hastings
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).

metropolis hastings eksempel independence sampler
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.
gibbs sampling
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)

gibbs sampling eks 4
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:
gibbs sampling kode
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

}

erfaringer med gibbs sampling
Erfaringer med Gibbs-sampling
  • Svært effektiv når parametrene er mer eller mindre uavhengige.
  • Mindre effektiv når det er stor avhengighet mellom parametrene.
  • Ikke alle parametre lar seg Gibbs-sample
  • Det går an å lage sampler av blokker av parametre ved å trekke noen via annen metodikk (independence sampler), deretter noen via Gibbs-aktig sampling og deretter ta Metropolis-Hastings-akseptans på hele blokken.
avledede parametre og prediksjoner
Avledede parametre og prediksjoner
  • Med ett sett trekninger av parametrene kan avledede parametre simpelthen beregnes og analyseres slik som parametrene selv.
  • Eks: prediksjon av vannføring: Q*~logN(,2). EQ*=exp(+0.52), der  er 2 hentet fra Gibbs-trekningene. Får EQ*473.7 m3/s.
modell valg
Modell-valg
  • Har et sett med modeller,Mk, og ønsker å beregne Pr(Mk|D)=f(D|Mk)Pr(Mk)/f(D).
  • Trenger derfor f(D|Mk). Modellene kan hver ha et parametersett k.
  • f(D|Mk) =  f(D| k , Mk) f(k |Mk) dk
  • Kan beregnes via Monte Carlo-metoder som importance sampling. Finn en forslagsfordeling gk(k) (f.eks. ved å tilpasse en multinormal fordeling til estimerte momenter fra MCMC-trekninger fra a’ posteriorifordelingen til modellparametrene) og beregn
  • Evt. finnes det en algoritme for å hoppe mellom modellrommene, kalt ”Reversible Jumps”, der andelen av iterasjoner der rutinen oppholder seg i modell k gir estimatet for Pr(Mk|D).
hierarkiske modeller
Hierarkiske modeller
  • Har så langt sett på tilfeller der parametrene lager data direkte. Men kan også ha parametre som gir andre parametre.
  • Vi har da flere lag av parametre og andre stokastiske variable (skjulte variable) før vi kommer til data. Bayes formel kan benyttes til å lære om lag k fra lag k-1.
hierarkiske modeller eks
Hierarkiske modeller - eks
  • Flomfrekvensanalyse: Antar at vi kan lære om hva fornuftige a’ priori-antagelser om GEV-parametrene for en stasjon j er ved å se på hele regionen:

hyper-parametre

 

 

 

1 1 1 GEV_parametre

2 2 2

3 3 3

Q1,1,Q1,2,…,Q1,n1

Q2,1,Q2,2,…,Q2,n2

Data

Q3,1,Q3,2,…,Q3,n1

forenklet bilde av hierarkisk modell
Forenklet bilde av hierarkisk modell

hyperparametre

     

j j j j=1,…,J=antall stasjoner

Qi,ji=1,…,nj Data

andre hierarkiske eksempler 1
Andre hierarkiske eksempler - 1
  • Vannføringskurve med dynamisk bunnvannstand: Qi=a(hi-h0,t)b, der h0,t= h0,t-1+t
  • Vannføringskurve der også vannstandsmålingene er å anse som uperfekte: Qi,real=a(hi,real-h0)b, Qi=Qi,real+ i, hi=hi,real+ i
  • Ifylling av hull med tidsseriemodell.
  • Vannførings-tidsserie med skjult regn-variabel: Qt=i aiRt-i + t, derRt= Rt-1+t
konklusjon
Konklusjon
  • Med Bayesiansk metodikk og muligheter for simulering, er det bare fantasien og programmeringsferdighetene som setter begrensningene!