Introduktion till Markovkedjor

(1)

Introduktion till

Markovkedjor

(2)

¨

Orebro universitet

Institutionen för naturvetenskap och teknik Matematik C, 76 – 90 högskolepoäng

Introduktion till Markovkedjor

Mattias Arvidsson

December 2011

Handledare: Holger Schellwat Examinator: Niklas Eriksen

Sj¨alvst¨andigt arbete, 15 hp Matematik, C–niv˚a, 76 – 90 hp

(3)

Sammanfattning

Uppsatsen är skriven med m˚alet att introducera läsaren till de grundläggande egenskaperna hos diskreta och ändliga Markovkedjor. Inneh˚allet är främst baserat p˚a källorna [6] och [7]. I kapitel 1, avsnitt 1.1 introduceras Markov-kedjor och deras grundläggande egenskaper och tolkningar, avsnitt 1.2 utvidgar de matematiska verktygen nödvändiga för att bevisa teorin kring Markov-kedjor och avsnitt 1.3 presenterar uppsatsens teoretiska tyngdpunkt genom satser och exempel. Uppsatsen avslutas med kapitel 2 där en tillämpning av Markovkedjor presenteras, kapitlet är mer informellt upplagt jämfört med kapitel 1 och syftar till att genom datorsimuleringar visa läsaren hur en Markovkedja kan användas.

(4)

(5)

Inneh˚

all

1 Markovkedjor 5

1.1 Introduktion till Markovkedjor . . . 5 1.2 En fixpunktssats f¨or affina avbildningar . . . 10 1.3 Irreducibla och regulj¨ara Markovkedjor . . . 13

2 Metropolis–Hastings algoritmen 19

2.1 Inledning . . . 19 2.2 M–H algoritmen för diskreta fördelningar . . . 21 2.3 M–H algoritmen för kontinuerliga fördelningar . . . 23

(6)

(7)

Kapitel 1

Markovkedjor

1.1 Introduktion till Markovkedjor

Vi har ett antal tillst˚and, S = (s0,s1,...,sr), d¨ar en process b¨orjar vid n˚agot

tillst˚and och successivt g˚ar fr˚an ett tillst˚and till nästa. Varje överg˚ang fr˚an ett tillst˚and till nästa kallas ett steg. Om kedjan är i tillst˚andet si g˚ar den

till tillst˚andet sj med sannolikhet pij, vilket kallas en ¨overg˚angssannolikhet.

Mängden av alla tillst˚and, utfallsrummet Ω, definieras av användaren. Om vi exempelvis är intresserade av vädret kan vi tänka oss att tillst˚anden si antar

utfallen regnigt (R), molnigt (M ) eller klart (K), vi f˚ar den diskreta utfalls-rummet si ∈ {R,M,K}. Om vi d˚a under en femdagarsperiod registrerar

R,R,R,M och K f˚ar vi S = (R,R,R,M,K).

Det b¨or n¨amnas att vi som standard definierar en vektor som en radvek-tor x = x1 · · · xn.

1.1.1 Definition (Markovegenskapen). En Markovkedja ¨ar en stokastisk process p˚a S med egenskapen

P (si+1|si) = P (si+1|si,si−1,...,s0) f¨or alla sk∈ S,

vilket inneb¨ar att sannolikheten pij enbart beror p˚a si. Den stokastiska

pro-cessen beror allts˚a enbart p˚a f¨oreg˚aende utfall.

Nedan ges definitionerna av stokastiska matriser och f¨ordelningsvektorer som anv¨ands genom uppsatsen.

1.1.2 Definition. En kvadratisk matris, P = [pij]n×n, d¨ar elementen i varje

radvektor ¨ar ¨overg˚angssannolikheter som kan summeras till ett, P_jpij =

1 för alla i, och 0 ≤ pij för alla i,j, kallas en stokastisk matris för en

Markovkedja. En radvektor, x, definierad som en radvektor i P , kallas en fördelningsvektor. En fördelningsvektor tillhör mängden K där

K= (

x= (x1· · · xn) ∈ Rn: xi ≥ 0 f¨or alla i och n X i=1 xi = 1 ) .

(8)

1.1.3 Definition. För en matris A = [aij] gäller A > 0 om aij > 0 för alla

i,j.

1.1.4 Lemma. F¨or tv˚a matriser A = [aij] och B = [bkl] s˚adana att A,B >

0, och matrisprodukten AB existerar g¨aller C := AB > 0.

Bevis. För varje element i C gällerP_rairbrl, där air > 0 och brl > 0 ger att

varje airbrl> 0. Vi f˚arPrairbrl > 0.

Nedan visas att produkten av tv˚a stokastiska matriser är en stokastisk matris, vi börjar med ett exempel för att illustrera den generella metoden. 1.1.5 Exempel. L˚at A och B vara stokastiska matriser enligt

A =   a11 a12 a13 a21 a22 a23 a31 a32 a33  , B =   b11 b12 b13 b21 b22 b23 b31 b32 b33  .

F¨orsta radvektorn av AB ¨ar

x₁= a₁₁b₁₁+ a₁₂b₂₁+ a₁₃b₃₁ a₁₁b₁₂+ a₁₂b₂₂+ a₁₃b₃₂ a₁₁b₁₃+ a₁₂b₂₃+ a₁₃b₃₃. Vi visar att x1 uppfyller definition 1.1.2:s krav p˚a en f¨ordelningsvektor.

0 ≤ aij ≤ 1 och 0 ≤ bij ≤ 1 för alla i,j medför 0 ≤ aijbij ≤ 1 för alla i,j. Vi

har nu kvar att visa att elementen i x1 kan summeras till ett.

a11b11+ a12b21+ a13b31+ a11b12+ a12b22+ a13b32+ a11b13+ a12b23+

a13b33 = a11(b11+ b12+ b13) + a12(b21+ b22+ b23) + a13(b31+ b32+ b33) =

a11· 1 + a12· 1 + a13· 1 = 1. Detta g¨aller f¨or alla radvektorer, xi, i matrisen

AB.

1.1.6 Sats. L˚at A = [aij]n×n och B = [bkl]n×n vara stokastiska matriser.

D˚a ¨ar matrisprodukten C := AB = [cij]n×n en stokastisk matris.

Bevis. 0 ≤ aij ≤ 1 och 0 ≤ bkl ≤ 1 medf¨or 0 ≤ aijbkl ≤ 1 f¨or alla i,j,k,l.

F¨or varje radvektor i matrisen C g¨aller ci = Prairbr1 . . . Prairbrn.

Summan av elementen i ci ¨ar P_rP_jairbrj =P_rairP_jbrj = 1 · 1 = 1 f¨or

alla i.

Genom sats 1.1.6 kan lemma 1.1.7 bevisas.

1.1.7 Lemma. L˚at P > 0 vara en stokastisk matris. D˚a ¨ar Pn _{> 0 en}

stokastisk matris f¨or alla n ∈ N+.

Bevis. Enligt sats 1.1.6 är P2 en stokastisk matris och enligt lemma 1.1.4 gäller P2 > 0. L˚at Pk > 0 för ett godtyckligt k ∈ N+ vara en stokastisk matris. D˚a är Pk_{P = P}k+1 _{en stokastisk matris s˚}_{adan att P}k+1_{> 0.}

Vi visar att produkten av en f¨ordelningsvektor och en stokastisk matris ¨

(9)

1.1.8 Sats. L˚at x vara en fördelningsvektor och l˚at P vara en stokastisk matris. D˚a gäller xP = x(1) där x(1) är en fördelningsvektor, dvs. för av-bildningen T : Rn→ Rn är x 7→ T (x) := xP linjär och T (x) ∈ K.

Bevis. L˚at x ∈ K och definiera avbildning T s˚adan att

T (x) = Pn_i=1pi1xi Pn_i=1pi2xi · · · Pn_i=1pinxi= xP.

T är en linjär transformation d˚a T (x + y) = (x + y)P = xP + yP = T (x) + T (y) och T (λx) = (λx)P = λ(xP ) = λT (x) för x,y ∈ K.

Vi visar att T (x) = x(1) ∈ K. 0 ≤ pij f¨or alla i,j och 0 ≤ xi f¨or

alla i ger 0 ≤ Pn_i=1pijxi f¨or alla j. Summan av alla j enheter i T (x) ¨ar

Pn j=1 Pn i=1pijxi = Pn i=1xi Pn j=1pij = 1 · 1 = 1.

Nedan följer ett första exempel p˚a en stokastisk matris. Det här exemplet ˚aterkommer genom större delen av uppsatsen.

1.1.9 Exempel. Överg˚angssannolikheterna utgör enheterna i en stokastisk matris P . Stokastiska matriser kan enkelt beskrivas genom ett exempel för väderprognoser. Säg att vädret för en dag kan vara regnigt (R), molnigt (M ) eller klart (K). Sannolikheten att vädret fr˚an en dag till en annan g˚ar fr˚an R till R är 1/2, fr˚an R till M är 1/4, fr˚an R till K är 1/4, fr˚an M till R är 1/2, fr˚an M till M är 0, fr˚an M till K är 1/2, fr˚an K till R är 1/4, fr˚an K till M är 1/4 och fr˚an K till K är 1/2. I matrisform blir det

P =   R M K R ,5 ,25 ,25 M ,5 0 ,5 K ,25 ,25 ,5  .

Sannolikheten att det exempelvis regnar tv˚a dagar fr˚an nu om det regnar i dag ¨ar

p(2)₁₁ = p(regnar om tv˚a dagar|valfritt v¨ader imorgon, regnar idag)

Det är skalär-produkten av matrisen P :s första rad och första kolumn. 1.1.10 Sats. L˚at P vara en stokastisk matris för en Markovkedja och l˚at n ∈ N+. Den ij:te enheten p(n)_ij i matrisen Pn ger sannolikheten att Markovked-jan, om den startar i tillst˚and si, är i tillst˚and sj efter n steg.

(10)

Bevis. L˚at P = [pij]n×n vara en stokastisk matris. Varje element, pij, ¨ar

sannolikheten att g˚a fr˚an ett tillst˚and till ett annat i ett steg, vilket skrivs som p(1)_ij . Enligt sats 1.1.6 ¨ar P2en stokastisk matris med radvektorer rioch

kolumnvektorer ci. Varje element i P2 ¨ar skal¨ar-produkten av matrisen P

med sig sj¨alv enligt

P2 =    r1· c1 · · · r1· cn .. . . .. ... rn· c1 · · · rn· cn    .

Vilket ¨ar sannolikheten, p(2)_ij , att g˚a fr˚an ett tillst˚and till ett annat i tv˚a steg. L˚at p(k)_ij vara sannolikheten att g˚a fr˚an tillst˚and si till tillst˚and sj i k

steg f¨or ett godtyckligt k ∈ N+. Sannolikheten att g˚a fr˚an ett tillst˚and till ett annat i tv˚a steg fr˚an Pk ¨ar p(k+1)_ij vilket ges av matrisen Pk+1. Genom induktion f˚ar vi att den ij:te enheten p(n)_ij i matrisen Pn_{ger sannolikheten att}

Markovkedjan, om den startar i tillst˚and si, ¨ar i tillst˚and sj efter n steg.

De sex första potenserna av den stokastiska matrisen för v˚ar väderprognos ¨

ar, med tre decimalers noggrannhet och utan ledande nollor, f¨oljande.

P1 ₌   R M K R ,500 ,250 ,250 M ,500 ,000 ,500 K ,250 ,250 ,500  , P2 ₌   R M K R ,438 ,188 ,375 M ,375 ,250 ,375 K ,375 ,188 ,438  , P3 ₌   R M K R ,406 ,203 ,391 M ,406 ,188 ,406 K ,391 ,203 ,406  , P4 ₌   R M K R ,402 ,199 ,398 M ,398 ,203 ,398 K ,398 ,199 ,402  , P5 =   R M K R ,400 ,200 ,399 M ,400 ,199 ,400 K ,399 ,200 ,400  , P6 =   R M K R ,400 ,200 ,400 M ,400 ,200 ,400 K ,400 ,200 ,400  .

För n större än sex är den stokastiska matrisen identisk med P6 _{för tre}

decimalers noggrannhet. Om det ¨ar regnigt idag (s0 = R) vad ¨ar d˚a

sanno-likheten att Markovkedjan ¨ar i tillst˚andet regnigt (s2 = R) om tv˚a dagar?

Enligt den andra potensen av den stokastiska matrisen, P2, är sannolikheten p(2)₁₁ = 0,438, vilket är sannolikheten att Markovkedjan är i tillst˚andet reg-nigt efter tv˚a steg om vi börjar i tillst˚andet regnigt (som tidigare beräknades med betingade sannolikheter).

1.1.11 Sats. L˚at P vara en stokastisk matris för en Markovkedja och l˚at xvara en fördelningsvektor som representerar startfördelningen. D˚a är san-nolikheten att Markovkedjan är i tillst˚and si efter n steg det i:te elementet i

(11)

vektorn

x(n)= xPn.

Bevis. Sannolikheten att Markovkedjan ¨ar i tillst˚and si efter ett steg ¨ar det

i:te elementet i vektorn xP . Enligt sats 1.1.10 ger Pnatt sannolikheten att Markovkedjan, om den startar i tillst˚and si, ¨ar i tillst˚and sj efter n steg. Vi

f˚ar att sannolikheten att en Markovkedja i tillst˚and si ¨ar i tillst˚and sj efter

n steg ¨ar det i:te elementet i vektorn x(n)= xPn.

Om vi som tidigare antar att vädret är regnigt idag kan x sättas till x= 1 0 0 Vi f˚ar x(1) = xP = 1 0 0   ,500 ,250 ,250 ,500 ,000 ,500 ,250 ,250 ,500  = ,5 ,25 ,25 x(2)= x(1)P = ,5 ,25 ,25   ,500 ,250 ,250 ,500 ,000 ,500 ,250 ,250 ,500  = ,438 ,188 ,375.

Vilket ¨ar ekvivalent med

x(2)= xP2 = 1 0 0   ,438 ,188 ,375 ,375 ,250 ,375 ,375 ,188 ,438  = ,438 ,188 ,375.

Som tidigare ¨ar sannolikheten att Markovkedjan ¨ar i tillst˚and regnigt efter tv˚a steg 0,438.

(12)

1.2 En fixpunktssats f¨

or affina avbildningar

Här ges den teori som behövs i avsnitt 1.3. Vi börjar med att definiera en fixpunkt.

1.2.1 Definition. En vektor x kallas en fixpunkt till en funktion T om T (x) = x.

1.2.2 Definition. Mängden K är konvex om för alla t ∈ [0,1] och för alla x,y ∈ K gäller det att tx + (1 − t)y ∈ K.

Genom definition 1.2.2 kan lemma 1.2.3 bevisas.

1.2.3 Lemma. L˚at n ∈ N+och l˚at mängden K vara en konvex mängd. För x₁_{, . . . ,x}_n_{∈ K gäller det att} 1

n(x1+ · · · + xn) ∈ K.

Bevis. Beviset utf¨ors genom induktion. L˚at k ∈ N+.

1. Induktionsantagande: x1, . . . ,xk∈ K ⇒ 1_k(x1+ · · · + xk) ∈ K

2. Induktionsp˚ast˚aende: x1, . . . ,xk+1∈ K ⇒ _k+11 (x1+ · · · + xk+1) ∈ K

Vi skriver om induktionsp˚ast˚aendet. y := _k+11 (x1+· · ·+xk+1) = _k+11 xk+1+ 1

k+1(x1+ · · · + xk) = k+11 xk+1+1+11 k

1

k(x1+ · · · + xk). Vi ska visa att y ∈ K,

s¨att z := 1_k(x1+ · · · + xk). Genom induktionsantagandet f˚as att z ligger

i K. Vi f˚ar y = _k+11 xk+1 + ₁₊11 k z. D˚a _k+11 ∈ [0,1] och 1 1+1 k ∈ [0,1] samt 1 k+1+1+11 k

= 1 om k ∈ N+ _f˚_{ar vi enligt definition 1.2.2 att y ∈ K.}

1.2.4 Definition. L˚at K vara en konvex mängd. En affin avbildning T : K→ K är en avbildning för x1,x2 ∈ K med egenskapen T (r1x1+ r2x2) =

r1T (x1) + r2T (x2) d¨arP2_i=1ri = 1 och ri ≥ 0.

Genom definition 1.2.4 kan sats 1.2.5 bevisas.

1.2.5 Sats. Om T : K → K är en affin avbildning följer det att för alla x₁_{, . . . ,x}_n_{∈ K och för alla r}₁_{, . . . ,r}_n _{∈ R där r}_i _{≥ 0 och} Pn

i=1ri = 1 g¨aller

det att T (Pn_i=1rixi) =Pn_i=1riT (xi) ∈ K f¨or n ∈ N+.

Bevis. Beviset utf¨ors genom induktion med utg˚angspunkt fr˚an beviset av sats 16.2 i [9]. L˚at k ∈ N+.

1. Induktionsantagande:Pk_i=1ri = 1 ⇒ T (Pk_i=1rixi) =Pk_i=1riT (xi)

2. Induktionsp˚ast˚aende: Pk+1_i=1 r_i′ = 1 ⇒ T (Pk+1_i=1r′_ixi) =Pk+1i=1 r

′

(13)

Beviset delas upp i tv˚a fall: r₁′ 6= 1 och r′₁= 1, för att undvika division med 0. Fallet d˚a r₁′ = 1 är trivialt, för fall r1′ 6= 1 f˚ar vi följande om vi bryter ut

r_k+1′ xk+1 T (r₁′x1+ . . . + r ′ k+1xk+1) = T       term1 z }| { r′_k+1xk+1+ (1 − r ′ k+1) term2 z }| { " r′₁ 1 − r′ k+1 x2+ . . . + r′_k 1 − r′ k+1 xk #       = r′_k+1T (xk+1) + (1 − r ′ k+1)T " r₁′ 1 − r′ k+1 x2+ . . . + r_k′ 1 − r′ k+1 xk #!

enligt definition 1.2.4. Vi kollar om P r′

i = 1 f¨or b˚ada termerna.

1. r′_k+1+ (1 − r′_k+1) = 1 2. r ′ 1 1−r_k+1′ + · · · + r′_k 1−r′_k+1 = r′1+···+r ′ k 1−r′_k+1 =hr′₁+ r′₂· · · + r′_k+1= 1 ⇔ r′₁· · · + r′_k= 1 − r′_k+1i= 1−r ′ k+1 1−r′_k+1 = 1

Av induktionsantagandet f¨oljer det att (1−r_k+1′ )T r′1 1−r_k+1′ x2+ . . . + r′_k 1−r′_k+1xk = (1 − r_k+1′ ) r1′ 1−r_k+1′ T (x2) + . . . + r′_k 1−r′_k+1T (xk) = r₁′T (x2) + · · · + r ′ kT (xk)

vilket sammanslaget med f¨orsta termen ger T (r′₁x1+r

′ 2x2+. . .+r ′ k+1xk+1) = r′ 1T (x2) + r ′ 2T (x2) + · · · + r ′ k+1T (xk+1).

1.2.6 Lemma. L˚at K vara en konvex mängd och l˚at T : K → K vara en linjär avbildning. D˚a är T en affin avbildning.

Bevis. S¨att x1, x2 ∈ K och r1, r2 ∈ R med r1+ r2 = 1. D˚a g¨aller

T (r1x1+ r2x2) = T (r1x1) + T (r2x2) = r1T (x1) + r2T (x2),

s˚a T ¨ar affin enligt definition 1.2.4 med till¨agget att ri≥ 0.

Följande sats presenteras utan bevis. Den intresserade läsaren hänvisas till [1].

1.2.7 Sats (Bolzano-Weierstrass). En begränsad talföljd i R inneh˚aller alltid en konvergent delföljd.

Genom sats 1.2.7 kan sats 1.2.8 bevisas. Beviset är avsett för vektorer i Rn. Notera även att en delmängd i Rn är kompakt om den är sluten och begränsad.

(14)

1.2.8 Sats (Markov-Kakutani). L˚at T vara en kontinuerlig affin avbildning av den icke-tomma, kompakta, konvexa m¨angden K till sig sj¨alv. D˚a har T en fixpunkt.

Bevis. L˚at z ∈ K och definiera T : K → K. Vi har T (z) ∈ K vilket ger z,T (z), T (T (z)), . . . , TN(z) ∈ K f¨or alla N ∈ N. Vi s¨oker en fixpunkt, x, s˚adan att {Tk_{(x) : k ∈ N} = {x} ⇒ T (x) = x. Definiera genomsnittet}

av N transformationer enligt xN = _N1 PN −1k=0 Tk(z). Enligt lemma 1.2.3 f˚ar

vi xN ∈ K. D˚a mängden K är kompakt är den sluten och en begränsad

delm¨angd i Rn_{. Enligt Bolzano-Weierstrass sats existerar det en konvergent}

delf¨oljd xNp av xN ∈ K som konvergerar mot ett x ∈ K s˚adana att

xNp− x

→ 0, p → ∞,

d˚a ¨ar x en fixpunkt. Det visas genom f¨oljande, m.h.a. sats 1.2.5.

xNp− T (xNp) = 1 Np Np−1 X k=0 Tk(z) − T   1 Np Np−1 X k=0 Tk(z)   = 1 Np Np−1 X k=0 Tk(z) − Np X k=1 Tk(z) = 1 Np z− TNp_(z) _{→ 0, p → ∞}

Där z− TNp_(z) _{≤ diam(K) är ändlig d˚}_{a K är begränsad. Det gäller}

¨ aven att T (xNp) − T (x) → 0, p → ∞ d˚a T ¨ar kontinuerlig. Vi f˚ar kx − T (x)k = limp→∞ xNp− T (xNp)

= 0 vilket ger att det existerar ett x s˚adant att T (x) = x.

(15)

1.3 Irreducibla och regulj¨

ara Markovkedjor

1.3.1 Definition. En Markovkedja kallas irreducibel om det ¨ar m¨ojligt att g˚a fr˚an varje tillst˚and till varje tillst˚and.

1.3.2 Definition. En Markovkedja kallas regulj¨ar om n˚agon potens av den stokastiska matrisen enbart har positiva element. Dvs. Pk > 0 f¨or n˚agot k ∈ N+.

Varje reguljär kedja är irreducibel men en irreducibel kedja behöver inte vara reguljär. I den stokastiska matrisen P för v˚ar väderprognos finns talet noll men exempelvis P2inneh˚aller enbart positiva element och fr˚an definition 1.3.2 är den stokastiska matrisen därför reguljär.

Skillnaden mellan irreducibla och regulj¨ara Markovkedjor kan visas genom f¨oljande exempel.

1.3.3 Exempel. L˚at en stokastisk matris vara s˚adan att P = 1 2 1 0 1 2 1 0 . F¨or udda potenser f˚ar vi P2n+1 =

0 1 1 0

och f¨or j¨amna potenser f˚ar vi P2n =

1 0 0 1

. Om vi börjar i tillst˚and 1 har vi Markovkedjan S = (1,2,1,2, . . .), det är allts˚a möjligt att g˚a fr˚an varje tillst˚and till varje tillst˚and s˚a Markovkedjan är irreducibel men ingen potens av P är s˚adan att Pk_{> 0}

s˚a Markovkedjan ¨ar inte regulj¨ar.

Sats 1.3.4 och 1.3.5 bevisas med utg˚angspunkt fr˚an [5] och [6]. Resultatet av sats 1.3.4 anv¨ands i beviset f¨or sats 1.3.5.

1.3.4 Sats. L˚at P > 0 vara en stokastisk matris. L˚at ε vara det minsta elementet i P . L˚at xT vara en kolumnvektor med st¨orsta element M0 och

minsta element m0. L˚at P xT vara en kolumnvektor med st¨orsta element M1

och minsta element m1. D˚a g¨aller M1≤ M0, m1 ≥ m0 och

M1− m1≤ (1 − 2ε)(M0− m0).

Bevis. L˚at yT _{vara vektorn x}T _{med alla element ersatta av M}

0 f¨orutom det

minsta m0 element. D˚a gäller yT ≥ xT. Varje element i P yT är som störst

om ε multipliceras med m0, dvs. det st¨orsta elementet, M1 i P yT m˚aste

vara begr¨ansat av

M1≤ εm0+ (1 − ε)M0.

Analogt erh˚alls en undre gräns för det minsta möjliga värdet (xT ersätts av yT med alla element ersatta av m0 förutom det största M0 elementet och ε

multipliceras med M0)

(16)

Vi f˚ar

M1− m1≤ εm0+ (1 − ε)M0− (εM0+ (1 − ε)m0) = (1 − 2ε)(M0− m0).

Genom sats 1.3.4 kan sats 1.3.5 bevisas.

1.3.5 Sats. L˚at P = [pij]n×n vara en regulj¨ar stokastisk matris f¨or en

Markovkedja. D˚a g¨aller

W = lim

n→∞P

n_,

där W är en matris med alla radvektorer w lika och positiva. Vi säger att P konvergerar mot matrisen W .

Bevis. En matris med alla radvektorer lika har identiska tal i varje kolumn. Vi bevisar existensen av W genom detta antagande. L˚at 0 < ε ≤ 1₂ ⇒ 0 ≤ (1 − 2ε) < 1 vara den stokastiska matrisens minsta element (d˚a P minst ¨ar av ordning 2 kan ε som st¨orst vara 1/2). L˚at xT _{vara en kolumnvektor med}

största värde M0 och minsta värde m0. L˚at Mn respektive mn vara

maxi-mum respektive minimaxi-mum elementen av vektorn Pn_xT_{. L˚}_{at M}

n−1respektive

mn−1 vara maximum respektive minimum elementen av vektorn Pn−1xT.

Fr˚an sats 1.3.4 f˚as

Mn− mn≤ (1 − 2ε)(Mn−1− mn−1)

f¨or n ≥ 1. Dessutom g¨aller det att M0 ≥ M1 ≥ . . . ≥ Mn och att m0 ≤

m1 ≤ . . . ≤ mn. Vi kan g¨ora omskrivningen

Mn− mn≤ (1 − 2ε)n(M0− m0) → 0, n → ∞.

Allts˚a g˚ar Pn_xT _{mot en kolumnvektor med alla element lika. Det ger att}

radvektorerna w i W ¨ar lika. Det kan visas genom att s¨atta xT _{= e}T_{, vi f˚}_ar

W eT = wT.

Vi utvidgar beviset till regulj¨ara matriser d¨ar ε′ _{≥ 0. Enligt definition}

1.3.2 existerar det ett Pk _s˚_{adant att ε}′ _{> 0, resten av beviset ¨ar analogt f¨or}

Pk, P2k, . . . , Pnk.

Nedan följer ett alternativt bevis för fallet d˚a kolumnvektorerna i Pk inte har alla element lika för n˚agot k utom i det fall d˚a P konvergerat mot W . Vi börjar med att bevisa ett lemma för viktade medelvärden. Varje fördelningsvektor uppfyller kraven p˚a vikterna i ett viktat medelvärde. 1.3.6 Lemma. L˚at ¯x = w1x1 + · · · + wnxn där 0 < wi < 1 för alla i,

P

iwi = 1 och det existerar minst ett xi s˚adant att xi 6= xj f¨or alla i 6= j.

(17)

Bevis. Vi f˚ar ¯x =Pn_i xiwi <Pni xmaxwi = xmaxPni wi= xmaxd˚a 0 < wi <

1 ochP_iwi = 1. Vi f˚ar maxixi > ¯x. Analogt f˚as att minixi< ¯x.

F¨or att beskriva metoden i det alternativa beviset b¨orjar vi med ett exempel.

1.3.7 Exempel. L˚at (P > 0) ∈ Mat (3,3) vara en stokastisk matris s˚adan att P =   p11 p12 p13 p21 p22 p23 p31 p32 p33  . Vi f˚ar P2 =   p11p11+ p12p21+ p13p31 p11p12+ p12p22+ p13p32 p11p13+ p12p23+ p13p33 p21p11+ p22p21+ p23p31 p21p12+ p22p22+ p23p32 p21p13+ p22p23+ p23p33 p31p11+ p32p21+ p33p31 p31p12+ p32p22+ p33p32 p31p13+ p32p23+ p33p33  .

Enligt lemma 1.3.6 f¨or en godtycklig kolumnvektor i P2 _{(v¨alj j godtyckligt)}

g¨aller f¨oljande.

minipij < p11p1j+ p12p2j+ p13p3j < maxipij

För ett godtyckligt element i P2 _{(välj i godtyckligt) gäller max}

i{pij} −

mini{pij} = ε1 > max{pi1p1j + pi2p2j + pi3p3j} − min{pi1p1j + pi2p2j +

pi3p3j} = ε2. Om vi forts¨atter resonemanget till P3, . . . , Pn f˚ar vi att ε1 >

ε2 > . . . > εn→ 0 d˚a n → ∞. Differensen mellan tv˚a element i en godtycklig

kolumn i Pn_g˚_{ar allts˚}_{a mot noll n¨ar n g˚}_{ar mot o¨andligheten. Vi f˚}_{ar matrisen}

W = lim n→∞P n₌   w1 w2 w3 w1 w2 w3 w1 w2 w3  , där w = w1 w2 w3 är en fördelningsvektor.

Nedan ges det alternativa beviset.

Bevis. Vi visar att en godtycklig kolumnvektor i W har alla element lika. L˚at P > 0 vara en stokastisk matris s˚adan att

P =    p11 · · · p1n .. . . .. ... pn1 · · · pnn    .

Enligt sats 1.1.6 ¨ar P2 en stokastisk matris. Vi f˚ar P2 =    Pn k=1p1kpk1 · · · Pnk=1p1kpkn .. . . .. ... Pn k=1pnkpk1 · · · Pn k=1pnkpkn    .

(18)

Enligt lemma 1.3.6 f¨or ett godtyckligt element i P2 _{g¨aller max}

i{pij} −

mini{pij} = ǫ1 > maxi{P_k=1pikpkj} − mini{P_k=1pikpkj} = ǫ2 f¨or alla

j. Analogt f˚as att ǫ1 > ǫ2 > . . . > ǫn = [maxi{pij} = mini{pij}] = 0 d˚a

n → ∞. Allts˚a g˚ar differensen mellan det st¨orsta och det minsta talet i varje kolumn i Pn mot noll.

I fallet d˚a P inneh˚aller n˚agot pij = 0 existerar det enligt definition 1.3.2

ett Pk_s˚_{adant att alla p}

ij > 0, resten av beviset ¨ar analogt f¨or Pk, P2k, . . . , Pnk.

I den stokastiska matrisen P för v˚ar väderprognos har vi att för tre decimalers noggrannhet W = lim n→∞P n₌   R M K R ,400 ,200 ,400 M ,400 ,200 ,400 K ,400 ,200 ,400  .

Sats 1.3.8 är en fundamental sats för reguljära Markovkedjor och tyngd-punkten av den här uppsatsen. I beviset används mycket av den teori vi upp till nu presenterat. Beviset utförs med utg˚angspunkt fr˚an [7] och [5]. 1.3.8 Sats. L˚at P vara en reguljär stokastisk matris, l˚at

W = lim

n→∞P

n

där w är den gemensamma radvektorn av W . D˚a gäller wP = w och en godtycklig radvektor v uppfyller vP = v om och endast om v är en multipel av w, dvs. vektorn w är unik för matrisen P och kallas en fixerad radvektor. Bevis. K är enligt definition 1.2.2 konvex. L˚at x ∈ K och definera avbild-ningen T : K → K enligt

T (x) = Pn_i=1pi1xi Pni=1pi2xi · · · Pni=1pinxi= xP

som enligt sats 1.1.8 ¨ar linj¨ar. Enligt sats 1.2.8 har T en fixpunkt w s˚adan att T (w) = w ⇔ wP = w.

L˚at v vara en godtycklig vektor s˚adan att vP = v. D˚a g¨aller vPn = v ⇒ vW = v. L˚at r vara summan av elementen i vektorn v, d˚a g¨aller vW = P_iviw1 · · · Piviwn = r w1 · · · wn = rw. Vi f˚ar v = rw,

allts˚a ¨ar v en multipel av w.

För den stokastiska matrisen P för v˚ar väderprognos har vi att

W =   R M K R ,400 ,200 ,400 M ,400 ,200 ,400 K ,400 ,200 ,400  .

(19)

Allts˚a ¨ar v˚ar fixerade radvektor f¨or P , w = ,4 ,2 ,4. Vi kan kon-trollera att w uppfyller villkoret i sats 1.3.8.

wP = ,4 ,2 ,4   ,5 ,25 ,25 ,5 0 ,5 ,25 ,25 ,5  = ,4 ,2 ,4

1.3.9 Sats. L˚at P vara en stokastisk matris för en reguljär Markovkedja med den fixerade radvektorn w. D˚a gäller det för en godtycklig fördelningsvektor x, xPn→ w d˚a n → ∞.

Bevis. Vi använder samma argumentation som i beviset för sats 1.3.8. L˚at x = x₁ · · · x_n vara en godtycklig fördelningsvektor. Enligt sats 1.3.8 gäller Pn→ W d˚a n → ∞, där W är en stokastisk matris s˚adan att

W =    w1 · · · wn .. . . .. ... w1 · · · wn    . Vi visar att xW = w. xW = P ixiw1 · · · Pixiwn = w1 · · · wn = w.

Vi kan allts˚a välja en godtycklig fördelningsvektor, säg x = ,2 0 ,8, och för W specificerad som innan f˚ar vi för v˚ar väderprognos

xW = ,2 0 ,8   ,400 ,200 ,400 ,400 ,200 ,400 ,400 ,200 ,400  = ,4 ,2 ,4.

Som ett sista exempel kan vi visa att w kan erh˚allas enligt sats 1.3.8.

w1 w2 w3   ,5 ,25 ,25 ,5 0 ,5 ,25 ,25 ,5  = w₁ w₂ w₃.

Vi f˚ar det linj¨ara ekvationssystemet

,5w1+ ,5w2+ ,25w3= w1

,25w1+ ,25w3 = w2

,25w1+ ,5w2+ ,5w3 = w3.

Sätt w1 = 1 (w1 kan väljas godtyckligt, det är enbart proportionerna mellan

w1, w2 och w3 som ¨ar intressanta). Vi f˚ar

,5 + ,5w2+ ,25w3= 1

,25 + ,25w3= w2

(20)

vilket har lösningen w = 1 ,5 1som ersatt med en proportionell fördelnings-vektor är w = ,4 ,2 ,4.

Slutligen bör det nämnas att w är vänster egenvektor med motsvarande egenvärde λ = 1 till den stokastiska matrisen P . Vi kan visa det genom att sätta wP = wλI ⇔ wP − wλI = 0 ⇔ w(P − λI) = 0.

För sats 1.3.9 utgörs egenrummet till egenvärde 1 av vektorerna som uppfyller {x : xW = x} som enligt beviset för sats 1.3.9 enbart är vektorn w, s˚a egenrummet är endimensionellt.

(21)

Kapitel 2

Metropolis–Hastings

algoritmen

2.1 Inledning

Vi bekantar oss nu med en tillämpning av Markovkedjor. Under 1900-talet har statistiker utvecklat metoder för att dra urval ur fördelningar som enbart ¨

ar möjliga att beräkna upp till en funktion proportionell mot den sökta. Flera algoritmer som approximerar en sökt fördelning har tagits fram; vi belyser Metropolis–Hastings (M–H) algoritmen.

Vi skiljer p˚a notationen för diskreta och kontinuerliga fördelningar. Teorin rörande M–H algoritmen är identisk för b˚ada fallen, skillnaden visar sig i deras tillämpningar. Vi börjar med att definiera en stokastisk variabel. 2.1.1 Definition. En stokastisk variabel, X, är en funktion som avbildar ett utfallsrum Ω p˚a R.

Till varje stokastisk variabel associerar vi i det diskreta fallet en sanno-likhetsfunktion och i det kontinuerliga fallet en t¨athetsfunktion.

2.1.2 Definition. En stokastisk variabel kallas diskret om den kan an-ta ett ¨andligt anan-tal olika v¨arden. Funktionen P (X = x) = pX(x) kallas

sannolikhetsfunktionen f¨or den stokastiska variabeln X om den uppfyller pX(x) ≥ 0 f¨or alla x ochP_xpX(x) = 1.

2.1.3 Definition. En stokastisk variabel kallas kontinuerlig om det exister-ar en funktion fX(x) s˚adan att

P (X ∈ A) = Z

A

fX(x)dx

för alla A ∈ R. Funktionen fX(x) kallas täthetsfunktionen för den stokastiska

(22)

I det diskreta fallet beskrevs sannolikheten att g˚a fr˚an ett tillst˚and x till nästa tillst˚and y av en stokastisk matris, P (x,y), av ändlig ordning med ett diskret utfallsrum. När vi talar om kontinuerliga funktioner ges överg˚ angs-sannolikheten fr˚an x till y av P (x,y) = P (Y ∈ A|X = x), som för varje x är en sannolikhetsfunktion dvs. P (x,y) ≥ 0 och RRP (x,y)dy = 1. Vi

betecknar en förslagsfördelning som den stokastiska matrisen Q(x,y) i det diskreta fallet och som täthetsfunktionen fY(y) i det kontinuerliga fallet. En

acceptanssannolikhetbetecknas som funktionen α(x,y).

M–H algoritmen beskrivs som ett matematisk verktyg och uppställnin-gen blir därför mer informell än tidigare, teorin som presenteras och deras tillämpningar är baserade p˚a [8], [3] och [4].

(23)

2.2 M–H algoritmen f¨

or diskreta f¨

ordelningar

Anta att vi vill simulera en sannolikhetsfunktion pX(x). V¨alj en regulj¨ar

stokastisk matris som förslagsfördelning, Q(x,y). M–H algoritmen ger en Markovkedja som konvergerar mot fördelningen av X enligt

P (x,y) = q(x,y)α(x,y) om x 6= y q(x,y) −Pu_6=xq(x,u)α(x,u) om x = y , f¨or acceptanssannolikheten α(x,y) = min 1,pX(y)q(y,x) pX(x)q(x,y) . Algoritmen kan beskrivas genom exempel 2.2.1.

2.2.1 Exempel.Vi vill simulera en sannolikhetsfunktion pX(x) = 1/6 2/6 3/6

med den stokastiska matrisen f¨or v˚ar v¨aderprognos i kapitel 1,

Q(x,y) =   R M K R 1/2 1/4 1/4 M 1/2 0 1/2 K 1/4 1/4 1/2  ,

som f¨orslagsf¨ordelning. Vi har tidigare definierat utfallsrummet som Ω = {R,M,K}. Vi f˚ar acceptanssannolikheten A(x,y) =      1 min1,pX(M )q(M,R) pX(R)q(R,M ) min1,pX(K)q(K,R) pX(R)q(M,K) min1, pX(R)q(R,M ) pX(M )q(M,R) 1 min1,pX(K)q(K,M ) pX(M )q(M,K) min1,pX(R)q(R,K) pX(K)q(K,R) min1,pX(M )q(M,K) pX(K)q(K,M ) 1      =      1 min1,2/6·1/2_1/6·1/4 min1,3/6·1/4_1/6·1/4 min1,1/6·1/4_2/6·1/2 1 min1,3/6·1/4_2/6·1/2 min1,1/6·1/4_3/6·1/4 min1,2/6·1/2_3/6·1/4 1      =   1 1 1 1/4 1 3/4 1/3 1 1  .

M–H algoritmen ger f¨oljande stokastiska matris P (x,y) =   1 − (1 · 1/4 + 1 · 1/4) 1 · 1/4 1 · 1/4 1/4 · 1/2 1 − (1/4 · 1/2 + 3/4 · 1/2) 3/4 · 1/2 1/3 · 1/4 1 · 1/4 1 − (1/3 · 1/4 + 1 · 1/4)   =   1/2 1/4 1/4 1/8 1/2 3/8 1/12 1/4 2/3  .

(24)

Vi kan enkelt se att pX(x)P = pX(x) och att lim n→∞P n₌   1/6 2/6 3/6 1/6 2/6 3/6 1/6 2/6 3/6  .

Fr˚an exempel 2.2.1 ser vi att det finns en risk f¨or division med noll f¨or vissa val av pX(x) och q(x,y).

(25)

2.3 M–H algoritmen f¨

or kontinuerliga f¨

ordelningar

Anta att vi vill simulera en täthetsfunktion fX(x). Välj en förslagsfördelning

fY(y) där värdemängden av Y är definitionsmängden av fX(x). M–H

algo-ritmen ger en Markovkedja som konvergerar mot f¨ordelningen av X enligt P (Yi,Yi+1) =

fY(Yi)α(Yi,Yi+1) om Yi6= Yi+1

fY(Yi) +R_Y_i+1_6=Y_ifY(Yi)(1 − α(Yi,Yi+1)) om Yi= Yi+1 ,

f¨or acceptanssannolikheten

α(Yi,Yi+1) = min

1,fX(Yi+1)fY(Yi) fX(Yi)fY(Yi+1)

.

Acceptanssannolikheten har för M–H algoritmen egenskapen att för en täthets-funktion som kan beräknas upp till en konstant, fX(x) = KgX(x) där

R

RgX(x)dx = 1 kommer K strykas i α(Yi,Yi+1) (vi s¨ager att t¨athetsfunktionen

inte beh¨over vara normerad).

Algoritmen kan beskrivas genom f¨oljande steg. 1. Vi b¨orjar i tillst˚and si = yi

2. Dra ett tal ur Y , vi f˚ar yi+1

3. Med sannolikhet α(yi, yi+1) ∈ [0,1] acceptera yi+1: si+1= yi+1

4. Annars stanna i nuvarande tillst˚and si+1= yi

5. Upprepa fr˚an steg 1 n antal g˚anger

2.3.1 Exempel. Vi kan illustrera algoritmen med följande exempel. Säg att vi vill simulera en X ∈ χ2₍₃₎–fördelning p˚a intervallet [0,10]. Vi f˚ar fX(x) = Kx3/2−1e−x/2 där K = 2

−3/2

Γ(3/2). L˚at den stokastiska variabeln Y

vara likformigt fördelad enligt Y ∈ U (0,10). Täthetsfunktionen för den lik-formiga fördelningen stryks i α(yi,yi+1). I simuleringsprogrammet R skriver

vi med utg˚angspunkt fr˚an [8] en algoritm som g¨or f¨oljande.

1. Dra en startpunkt ur Y . Vi f˚ar y0 = 7,627578, vilket ger S = (7,627578)

2. Dra ett tal ur Y , vi f˚ar y1 = 1,151297

3. Med sannolikhet α(y0,y1) = min

1,fX(y1)fY(y0) fX(y0)fY(y1) acceptera 1,151297, vi f˚ar:

(26)

Figur 2.1: Exempel p˚a stegen i en Markovkedja genererad av en M-H algo-ritm, referens [2]. α(7,627578; 1,151297) = min 1,K1,151297 3/2−1_e−1,151297/2 K7,6275783/2−1_e−7,627578/2 ! = min 1,1,151297 3/2−1_e−1,151297/2 7,6275783/2−1_e−7,627578/2 ! = min(1; 9,901635) = 1

4. Vi accepterar 1,151297 med 100% sannolikhet, vilket ger S = (7,627578; 1,151297) 5. Upprepa proceduren n antal g˚anger. Resultatet blir att S konvergerar

mot f¨ordelningen av X d˚a n → ∞

Resultatet av myMCMC2(102) ges i figur 2.2. Vid ett st¨orre antal itera-tioner ger S, genom P (x,y), dragningar ur en χ2₍₃₎–f¨ordelning. Programkoden redovisas nedan.

(27)

Histogram of X[−seq(nsim * 0.25, nsim)] X[−seq(nsim * 0.25, nsim)] Frequency 0 2 4 6 8 10 0 10 20 30 40 50 60

Histogram of X[−seq(nsim * 0.5, nsim)]

X[−seq(nsim * 0.5, nsim)] Frequency 0 2 4 6 8 10 0 20 40 60 80 100

Histogram of X[−seq(nsim * 0.75, nsim)]

X[−seq(nsim * 0.75, nsim)] Frequency 0 2 4 6 8 10 0 50 100 150 Histogram of X X Frequency 0 2 4 6 8 10 0 50 100 150 200

Figur 2.2: Histogram av simuleringen för 102· 0.25, 102· 0.5, 102· 0.75 och 102 iterationer, röd linje är den sanna fördelningen.

myMCMC2 <- function(nsim = 10^4){

X <- rep(runif(n = 1, min = 0, max = 1),nsim) # initialize the chain Z <- dchisq(seq(0, 10, length.out = nsim), df = 3)

for (i in 2:nsim){

Y <- runif(n = 1, min = 0, max = 10)

rho <- min(1,dchisq(Y, df = 3)*dunif(X[i-1], min = 0, max = 10)/ (dchisq(X[i-1], df = 3)*dunif(Y, min = 0, max = 10)))

X[i] <- X[i-1] + (Y-X[i-1])*(runif(n = 1, min = 0, max = 1)<rho) }

attach(mtcars) par(mfrow=c(2,2))

hist(X[-seq(nsim*.25, nsim)]) par(new=TRUE)

plot(Z, yaxt=’n’,xaxt=’n’, ann=FALSE,col="red") hist(X[-seq(nsim*.5, nsim)])

par(new=TRUE)

(28)

hist(X[-seq(nsim*.75, nsim)]) par(new=TRUE)

plot(Z, yaxt=’n’,xaxt=’n’, ann=FALSE,col="red") hist(X)

par(new=TRUE)

plot(Z, yaxt=’n’,xaxt=’n’, ann=FALSE,col="red") return()

(29)

Litteraturf¨

orteckning

[1] D.A. Brannan. A first course in mathematical analysis. Cambridge Univ Pr, 2006.

[2] I. Cosma. Markov chains and monte carlo methods. A.I.M.S, 2010. [3] P. Diaconis. The markov chain monte carlo revolution. Bull. Amer.

Math. Soc.(NS), 46(2):179–205, 2009.

[4] R. Grey. Advanced statistical computing. Department of Biostatistics and Computational Biology, Harvard, 2002.

[5] C.M. Grinstead and J.L. Snell. Introduction to probability. Amer Math-ematical Society, 1997.

[6] J.G. Kemeny and J.L. Snell. Finite markov chains. Springer, 1960. [7] G. Maltese. A simple proof of the fundamental theorem of finite markov

chains. The American mathematical monthly, 93(8):629–630, 1986. [8] C.P. Robert and G. Casella. Introducing Monte Carlo Methods with R.

Springer Verlag, 2010.