Aanvullend dictaat Stochastische Operations Research I. H.C. Tijms

Aanvullend dictaat Stochastische Operations Research I H.C. Tijms november 2008

2

Contents 1 Markov Keten Monte Carlo Methoden 1.0.1 Reversibele Markov ketens . . . 1.0.2 Metropolis-Hastings algoritme . 1.0.3 De Gibbs sampler . . . . . . . . 1.1 Opgaven . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

2 Stochastic Dynamic Programming 2.1 Stochastische dynamische programmering . . 2.1.1 Een dobbelspel en optimaal stoppen 2.1.2 Het spel rood en zwart . . . . . . . . 2.2 Investeringsprobleem en de Kelly strategie . 2.2.1 Dynamische programmering met een ningsduur . . . . . . . . . . . . . . . 2.3 Opgaven . . . . . . . . . . . . . . . . . . . . 3 4

Appendix 1: Poisson process

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . stochastische . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . . . . . . . . . . plan. . . . . . . .

7 7 10 16 19 21 21 21 24 25 28 30 33

Appendix 2: Renewal-reward processes 43 4.1 Basic theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2 Poisson arrivals see time averages . . . . . . . . . . . . . . . . . . 48 4.3 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3

4

CONTENTS

Voorwoord Voor u ligt een aanvullend diktaat voor het college Stochastische Operations Research I. Het onderwerp discrete-tijds Markov ketens staat centraal in het college SOR I. Dit onderwerp wordt op een meer elementair niveau behandeld in hoofdstuk 15 van het eerder bij het tweedejaars college Simulatie gebruikte boek H.C. Tijms, Understanding Probability en op een hoger niveau in het boek H.C. Tijms, A First Course to Stochastic Models. Dit laatstgenoemde boek zal het te gebruiken boek zijn bij het vervolgcollege Stochastische Operations Research II waarin onder meer continue-tijds Markov ketens met toepassingen in de wachttijdtheorie aan de orde komen. In het college SOR I zullen we ons echter voor de discrete-tijds Markov ketens baseren op het materiaal in hoofdstuk 15 van het boek Understanding Probability en dit aanvullen met materiaal over Markov keten Monte Carlo methoden, een methode die steeds meer in allerlei gebieden gebruikt wordt (hoofdstuk 1 uit dit aanvullend diktaat). Verder wordt in het college SOR I een uit het boek Operationele analyse aangepast hoofdstuk over stochastische dynamische programmering (in feite Markov ketens met besturing) behandeld en dit is hoofdstuk 2 in het aanvullende diktaat. In het college dat zes weken lang twee keer twee uur college per week beslaat, zal sprake zijn van een mengvorm van hoor- en werkcollege. De opzet zal zo zijn dat iedere deelnemer actief betrokken wordt bij het maken van de opgaven. De ervaring wijst uit dat tentamen normaal gesproken niet haalbaar is indien niet actief deelgenomen aan de werkvorm. Verder is op vrijwillige en individuele basis een bonus van maximaal 0.5 punt te verkrijgen, zie hieronder bij college 4. College 1 (a) Het Markov-keten model, paragraaf 15.1 (uit boek Understanding Probability, 2de druk, Cambridge University Press, 2007). (b) Opgaven 15.1, 15.2 15.4, 15.5 en 15.6. College 2 (a) Tijdsafhankelijk gedrag van Markov-ketens en absorberende Markovketens, de paragrafen 15.2 en 15.3. (b) Opgaven 15.11, 15.12, 15.14, 15.18, 15.20 en 15.24. College 3 (a) Evenwichtsanalyse van Markov-ketens, paragraaf 15.4. (b) Opgaven 15.25, 15.29, 15.32 en 15.34 College 4 Hoofdstuk 1 over Markov keten Monte Carlo Methoden uit aanvullend diktaat.

6

CONTENTS

NB: Aan de opgaven 2 en 3 kan op vrijwillige basis gewerkt worden en dit geeft dan een bonus van maximaal 0.5 punt (inleveren dient te geschieden uiterlijk op het laatste college; lengte verslag tussen 5 en 10 bladzijden exclusief de computerprogramma’s die als bijlage toegevoegd dienen te worden). College 5 (a) Stochastische dynamische programmering, paragraaf 2.1 uit aanvullend diktaat. (b) Opgaven 1, 2, 3 en 4. College 6 (a) Stochastische Dynamische programmering, paragraaf 2.2 uit aanvullend diktaat. (b) Opgaven 6, 7, 8 en 9.

Chapter 1 Markov Keten Monte Carlo Methoden De Markov keten Monte Carlo methode is een methode om uit een (multivariate) kansdichtheid π(x) te simuleren die op een multiplicative constante na bekend is, waarbij het niet mogelijk is de multiplicatieve constante rechtstreeks door normalisatie te berekenen. Deze krachtige methode wordt in de praktijk veelvuldig gebruikt, met name in de Bayesiaanse statistiek die met deze simulatiemethode veel aan toepasbaarheid heeft gewonnen. Het basisidee van de simulatiemethode is simpel. In feite komt de methode neer op het construeren van een Markov keten waaruit direct random trekkingen gedaan kunnen worden en die π(x) als evenwichtsverdeling heeft. De theoretische basis van zo’n constructie is het begrip reversibele Markov keten dat we eerst zullen bespreken.

1.0.1

Reversibele Markov ketens

Het begrip reversibiliteit introduceren we aan de hand van het Ehrenfest model. In dit model is sprake van twee compartimenten A en B die tezamen r deeltjes bevatten. Elke keer wordt één van de deeltjes random gekozen en van compartiment verwisseld. Definiëren we Xn als het aantal deeltjes in compartiment A na de nde stap, dan is {Xn } een Markov keten met toestandsruimte I = {0, 1, . . . , r} en 1-staps overgangskansen pi,i−1 = ri , pi,i+1 = r−i en pij = 0 anders. De r evenwichtsvergelijkingen πj = pj−1,j πj−1 + pj+1,j πj+1 , waarbij π0 = p10 π1 en πr = pr−1,r πr−1 , kunnen worden beschreven tot (ga na!): pj,j−1 πj = pj−1,j πj−1

voor j = 1, . . . , r.

Bedenken we dat pij = 0 voor |j −i| > 1, dan volgt uit deze relatie dat de Markov keten {Xn } de eigenschap heeft dat pjk πj = pkj πk voor alle j, k ∈ I. Definitie 7.5. Een irreducibele Markov keten 1

1

{Xn } met een eindige toestand-

Een Markov keten heet irreducibel als elk tweetal toestanden onderling bereikbaar zijn, (n) d.w.z. voor elke i en j is er een n met pij > 0.

7

8

CHAPTER 1. MARKOV KETEN MONTE CARLO METHODEN

sruimte I heet reversibel als voor de unieke evenwichtsverdeling {πj } van de Markov keten geldt dat πj pjk = πk pkj

voor alle j, k ∈ I.

In woorden, in de evenwichtsituatie is het verloop van het proces terugwaarts in de tijd probabilistisch gezien identiek aan het verloop van het proces voorwaarts in de tijd. Het bewijs is simpel. Voor elke n ≥ 1 geldt P (Xn = j) = πj voor alle j, zoals bij definitie 7.2 beargumenteerd is. Dus P (Xn−1 = j | Xn = k) =

P (Xn = k | Xn−1 = j)P (Xn−1 = j) pjk πj . = P (Xn = k) πk

Dit geeft P (Xn−1 = j | Xn = k) = pkj = P (Xn = j | Xn−1 = k), omdat pjk πj = pkj op grond van reversibiliteit. πk In het vervolg zullen we ons beperken tot irreducibele Markov ketens met een eindige toestandsruimte. Een belangrijk resultaat dat we herhaaldelijk zullen gebruiken, is het volgende: Stelling 7.5 Laat {Xn } een irreducibele Markov keten zijn met een eindige toestandsruimte I. Als er een kansverdeling {aj , j ∈ I} is zodat aj pjk = ak pkj

voor alle j, k ∈ I,

dan geldt aj = πj voor alle j ∈ I, waarbij {πj } de unieke evenwichtsverdeling van de Markov keten is. Bewijs. Sommeren we de vergelijking voor ajP over k ∈ I, dan vinden we met P k∈I pjk = 1 de evenwichtsvergelijkingen aj = k∈I ak pkj voor j ∈ I. Op grond van stelling 7.3 is de evenwichtsverdeling {πj } de enige kansverdeling die hieraan voldoet, waarmee de stelling bewezen is. Een interessante vraag is de volgende. Stel dat {aj , j ∈ I} een kansmassa functie op een eindige verzameling I is met aj > 0 voor alle j. Is het mogelijk een Markov keten te construren die {aj , j ∈ I} als unieke evenwichtsverdeling heeft? Het antwoord is ja. In feite, zijn zelfs oneindig veel Markov ketens met deze evenwichtsverdeling te construeren. De fysische constructie van zo’n Markov keten is als volgt: stel dat de huidige toestand j is, kies dan random één van de andere toestanden, zeg toestand k. Deze toestand k is de volgende toestand van de Markov keten als ak > aj , anders is toestand k met kans ak /aj de volgende toestand van de Markov keten en blijft de Markov keten in de huidige toestand j met kans 1 − ak /aj . Dus, met N = |I| het aantal toestanden in I, definieer een Markov keten op I met de 1-staps overgangskansen: voor elke j, k ∈ I met j 6= k geldt aj pjk =

1 1 min(ak , aj ) = ak min(aj /ak , 1) = ak pkj . N −1 N −1

9 Deze Markov keten is irreducibel en aperiodiek (ga na!). Verder voldoet de Markov keten aan de reversibiliteitsconditie. Stelling 7.5 geeft nu dat de Markov keten {aj } als unieke evenwichtsverdeling heeft. Opmerking: {aj } is ook de unieke evenwichtsverdeling is van de Markov keten die we verkrijgen door in de bovenstaande pjk ’s de constante N 1−1 te vervangen door een constante γ > 0 met (N − 1)γ ≤ 1.

Simulated annealing algoritme Reversibele Markov ketens liggen ten grondslag aan het simulated annealing algoritme. Dit is een zoekmethode om het absolute minimum van een veelal gecompliceerde functie te bepalen op een eindig maar zeer groot domein. Het kernidee van het algoritme is om volgens een kansverdeling van het ene punt naar het andere punt te bewegen zodat de zoekprocedure ook uit een locaal minimum kan ontsnappen. Wij geven slechts het ruwe idee van het algoritme. Stel dat c(i) een gegeven functie is op een eindige verzameling I. Voor elk punt i ∈ I is een locale omgevingsverzameling N (i) van punten gekozen met i ∈ / N (i) zodanig dat j ∈ N (k) als k ∈ N (j). Verder wordt verondersteld dat voor elk tweetal verschillende punten l en k er een keten l0 = l, l1 , . . . , lr = k is met lv ∈ N (lv−1 ) voor v = 1, . . . , r. Voor het gemak nemen we aan dat elke N (j) uit eenzelfde aantal punten bestaat. Wij definiëren nu als volgt een Markov keten op I. Als de huidige toestand van de Markov keten j is, dan wordt random een kandidaattoestand k uit N (j) gekozen. De volgende toestand van de Markov keten is gelijk aan k als c(k) < c(j); anders is de volgende toestand gelijk aan k met kans e−c(k)/T /e−c(j)/T en gelijk aan de huidige toestand j met kans 1−e−c(k)/T /e−c(j)/T . Hierbij is T > 0 een besturingsparameter. Met andere woorden, de Markov keten is gedefinieerd door de 1-stapsovergangskansen  1  M min(e−c(k)/T /e−c(j)/T , 1) voor k ∈ N (j) pjk = P  1 − l6=j pjl voor k = j en pjk = 0 anders, waarbij M is het aantal elementen in elk van de N (j)’s. Voor deze irreducibele Markov keten geldt e−c(j)/T pjk = e−c(k)/T pkj

∗

voor alle j, k ∈ I,

oftewel de Markov keten is reversibel. Het bewijs is simpel. Neem j. Voor ´ ³ k 6= e−c(k)/T −c(j)/T 1 k∈ / N (j) geldt pjk = pkj = 0 en voor k ∈ N (j) is e = min 1, M e−c(j)/T ³ −c(j)/T ´ 1 min(e−c(j)/T , e−c(k)/T ) = e−c(k)/T M1 min 1, ee−c(k)/T . De Markov keten heeft M dus de evenwichtskansen X 1 πi = e−c(i)/T voor i ∈ I met A = e−c(k)/T . A k∈I ∗∗

Als |N (j)| niet hetzelfde is voor alle j, neem dan M = maxj |N (j)|.

10


Als de functie c(i) het P absoluut minimum aanneemt in een uniek punt m, dan volgt uit πm = 1/(1 + Pk6=m e−(c(k)−c(m))/T ) dat πm → 1 als T → 0 (ga zelf na dat algemeen geldt dat i∈M πi → 1 als T → 0 met M is de verzameling van de punten waarin de functie c(i) het absolute minimum aanneemt). Dit resultaat is de basis van het simulated annealing algoritme. In dit algoritme wordt met een grotere waarde voor de zogenoemde afkoeltemperatuur T begonnen en laat men in elke iteratiestap n van het algoritme de parameter T dalen, bijvoorbeeld volgens T =C/ln(n + 1) met C > 0 een constante. Het algoritme wordt gestart met een toestand X0 = i0 . Achtereenvolgens worden dan toestanden X1 = i1 , X2 = i2 , . . . , XN = iN gegenereerd met N groot en het absolute minimum van de functie c(i) wordt dan geschat door mink=0,1,...,N c(ik ). Een interessante toepassing van het simulated annealing algoritme is het handelsreizigersprobleem. Stel dat uitgaande van stad 0 de steden 1, . . . , r bezocht moeten worden, waarbij elk van deze steden slechts één keer aangedaan mag worden en weer teruggekeerd moet worden in de beginstad 0. Veronderstel dat niet-negatieve kosten c(i, j) worden gemaakt wanneer vanuit stad i de stad j als volgende stad wordt aangedaan. Een permutatie x = (x1 , . . . , xr ) van de gehele getallen 1, . . . , r geeft dan een route met de interpretatie dat vanuit stad 0 naar stad x1 wordt gegaan, vanuit stad x1 naar stad x2 , etc en uiteindelijk Pr+1vanuit stad xr weer terug naar stad 0. De kosten van zo’n route zijn c(x) = i=1 c(xi−1 , xi ) met x0 = xr+1 = 0. In het simulated annealing algoritme zou je voor elke route x als omgevingsverzameling N (x) kunnen kiezen al die routes die ontstaan door de verwisseling van twee elementen in de permutatie (x1 , . .¡. ¢, xr ). Iedere omgevingsverzameling N (x) bestaat dan uit hetzelfde aantal van 2r elementen.

1.0.2

Metropolis-Hastings algoritme

In Markov keten theorie is het gebruikelijk een evenwichtsverdeling te zoeken voor een Markov keten. Bij de Markov keten Monte Carlo methode is in feite het omgekeerde het geval: op een multiplicatieve constante na, is een kansverdeling π(s) op een eindige doch zeer grote waardenverzameling S gegeven met π(s) > 0 voor alle s. Het doel is om een Markov {Xn } keten te vinden die de kansverdeling {π(s)} als unieke evenwichtsverdeling heeft en waaruit we met simulatie een P prestatiemaat van het type s∈S f (s)π(s) voor een gegeven functie f op S kunnen schatten, gebruikmakend van de ergodenstelling n

X 1X lim f (Xk ) = f (s)π(s) n→∞ n s∈S k=1

met kans 1.

Hoe de 1-stapsovergangskansen van zo’n Markov keten te vinden? Noteer deze overgangskansen als pM H (s, t). Het idee is om eerst zogenoemde kandidaatovergangskansen te kiezen en die vervolgens zodanig aan te passen dat de reversibiliteitsconditie voor Markov ketens vervuld is. Voor elke s ∈ S, kies je

11 een kansmassa functie (discrete kansdichtheid) q(t | s), t ∈ S. Dit doe je zodanig dat de Markov keten met q(t | s) als 1-staps overgangskans van toestand s naar toestand t irreducibel is. Ingeval de kandidaat-kansen q(t | s) meteen al zouden voldoen aan π(s)q(t | s) = π(t)q(s | t) voor alle s, t ∈ S. Dan is de reversibiliteitsconditie vervuld en kunnen we op grond van stelling 7.5 stellen dat {π(s)} de evenwichtsverdeling is van de Markov keten met de q(t | s) als 1-staps overgangskansen en hebben we de gezochte Markov keten geconstrueerd. In het algemeen zal niet voor alle s, t de bovenstaande gelijkheid gelden. Stel dat voor de combinatie (s, t) met s 6= t het ongelijkheidsteken geldt. Neem zonder beperking aan dat π(s)q(t | s) > π(t)q(s | t). In dit geval kunnen we, losjes gezegd, stellen dat het proces te frequent van s naar t gaat en te weinig vaak van t naar s. Een geschikte manier om dit te herstellen is door het aantal transities van s naar t te verminderen door middel van een zogenoemde acceptatiekans α(s, t): met kans 1 − α(s, t) vindt de transitie van s naar t geen doorgang en blijft het proces in toestand s. De gezochte 1-staps overgangskansen pM H (s, t) worden dan gekozen als: pM H (s, t) = q(t | s)α(s, t) voor s 6= t. Hoe α(s, t) en α(t, s) te kiezen? De bovenstaande ongelijkheid vertelt ons dat transities van t naar s niet vaak genoeg gebeuren zodat het logisch is om in elk geval α(t, s) = 1 te kiezen. De keuze van α(s, t) wordt bepaald door de wens dat de 1-staps overgangskansen pM H (s, t) = q(t | s)α(s, t) voldoen aan de reversibiliteitsconditie. Dit geeft met α(t, s) = 1 de eis π(s)q(t | s)α(s, t) = π(t)q(s | t)α(t, s) = π(t)q(s | t). Oftewel, met α(s, t) =

π(t)q(s | t) π(s)q(t | s)

kunnen we bovenstaande ongelijkheid aanpassen tot een gelijkheid. Samenvattend, bij kandidaat-kansen q(t | s) kiezen we de acceptatiekansen α(s, t) volgens · ¸ π(t)q(s | t) α(s, t) = min ,1 voor alle s, t ∈ I. π(s)q(t | s) De 1-staps overgangskansen pM H (s, t) van de gezochte Markov keten definiëren we door

½ pM H (s, t) =

q(t | P s)α(s, t) voor t 6= s 1 − t6=s q(t | s)α(s, t) voor t = s.

12


De Markov keten met deze 1-staps overgangskansen voldoet aan de reversibiliteitsconditie π(s)pM H (s, t) = π(t)pM H (t, s)

voor alle s, t.

De aanname is gemaakt dat de kandidaat Markov keten met de 1-staps overgangskansen q(t | s) irreducibel is. Deze aanname is essentieel en impliceert dat de Markov keten met de 1-staps overgangskansen pM H (s, t) ook irreducibel is zodat de kansverdeling π(s) de unieke evenwichtsverdeling is van deze Markov keten.∗ Een belangrijke opmerking is dat voor de constructie van deze gezochte Markov keten het voldoende is de kansen π(s) te kennen op een multiplicatieve constante na. Immers in α(s, t) wordt alleen het quotiënt π(t)/π(s) gebruikt! Het toestandsverloop in de Markov keten met de pM H (s, t) als 1-staps overgangskansen is als volgt te beschrijven. Als de huidige toestand s0 is, dan wordt een kandidaat-toestand t1 geloot volgens de kansdichtheid {q(t | s0 ), t ∈ S}. Deze kandidaat-toestand wordt met kans α(s0 , t1 ) geaccepteerd als volgende toestand s1 van de Markov keten, terwijl met kans 1 − α(s0 , t1 ) de kandidaat-toestand t1 wordt verworpen in welk geval de volgende toestand s1 van de Markov keten gelijk is aan de huidige toestand s0 .

Metropolis-Hastings algoritme Dit algoritme genereert een rij opeenvolgende toestanden uit een Markov keten {Xn } die {π(s), s ∈ S} als unieke evenwichtsverdeling heeft, waarbij de kansen π(s) > 0 op een multiplicatieve constante na gegeven zijn. Stap 0. Kies voor elke s ∈ S een kansdichtheid q(t | s). Neem een begintoestand s0 uit S. Laat X0 := s0 en n := 1. Stap 1. Trek een kandidaat-toestand tn uit de kansdichtheid {q(t | sn−1 ), t ∈ S}. Bereken de acceptatiekans · ¸ π(tn )q(sn−1 |tn ) ,1 . α = min π(sn−1 )q(tn | sn−1 ) Stap 2. Trek een random getal u uit (0, 1). Als u ≤ α, dan accepteer tn en sn := tn ; anders, sn := sn−1 . Stap 3. Xn := sn en n := n + 1. Herhaal stap 1 met sn−1 vervangen door sn . ∗∗

Als voor de Markov keten met de q(t | s) als 1-stapsovergangskansen de aanname van irreducibiliteit was verzwakt tot geen twee disjuncte fuiken, dan zou de Markov keten met de 1-staps overgangskansen (pM H (s, t)) meerdere disjuncte fuiken kunnen hebben. Dit is het geval als voor een zekere ta geldt q(ta | s) = 1 voor elke s; de acceptatiekans α(s, ta ) is dan nul voor elke s 6= ta zodat elke s 6= ta een absorberende toestand is.

13 Als de gekozen kansdichtheden q(t | s) symmetrisch zijn (d.w.z., q(t | s) = q(s | t) voor³alle s, t), ´dan vereenvoudigt in het algoritme de acceptatiekans α tot α = π(tn ) min π(s , 1 , hetgeen het oorspronkelijke algoritme van Metropolis was. n−1 ) P Stel dat je voor gegeven functie f de getalwaarde E[f (X)] = s∈S f (s)π(s) wilt berekenen met X een stochast die π(s) als kansdichtheid heeft. Pas dan m stappen van een Metropolis-Hastings algoritme toe met m voldoende groot. Op grond van de ergodenstelling voor Markov ketens schat je E[f (X)] met de gegenereerde rij van toestanden s1 , s2 , . . . , sm door m

1 X f (sk ). m k=1 Het Metropolis-Hastings algoritme is beschreven voor de situatie met een discrete toestandsverzameling, waarbij je dus voor toestand s de kandidaat-toestand t trekt uit de discrete kansdichtheid q(t | s), t ∈ S. Een nadere bestudering van het algoritme leert dat het algoritme woordelijk doorgaat voor de situatie met een continue toestandsverzameling mits enkele netheidscondities vervuld zijn. In die situatie trek je elke keer een kandidaat-toestand t uit een continue kansdichtheid. Bij de implementatie van het Metropolis-Hastings algoritme is de vraag hoe de kandidaat-dichtheden q(t | s) te kiezen. Vele keuzes zijn mogelijk en als de ene keuze niet bevredigend werkt dan kan een andere keuze geprobeerd worden. De keuze van de begintoestand kan ook van invloed zijn. Het implementatie aspect wordt uitgebreid besproken in het klassieke artikel S. Chib and E. Greenberg, Understanding the Metropolis-Hastings algorithm, The American Statistician, Vol 49 (1995), blz 327-335. Een goede menging in de Markov-keten is essentieel, d.w.z. de Markov keten beweegt voldoende snel door het gehele waardenbereik van de te simuleren dichtheid π(x). Ongewenst is dus dat de Markov keten lange tijd in dezelfde toestand blijft zoals bij een lage acceptatiekans het geval is. Veelgebruikte keuzes zijn: (a) de onafhankelijkheidskeuze waarin q(t | s) gegeven wordt door een kansdichtheid q(t) die onafhankelijk van s is (het is belangrijk dat de staart van de kansdichtheid q(x) die van de te simuleren kansdichtheid π(x) domineert). (b) de random-walk keuze waarin de nieuwe kandidaat-toestand t bepaald wordt door t = s+Z met Z een stochast die een gegeven kansverdeling heeft (in dit geval hangt de kansdichtheidsfunctie q(t | s) alleen af van t − s). Als de kansdichtheid van Z symmetrisch is rond de oorsprong, dan geldt q(t | s) = q(s | t) voor alle s, t zodat in het algoritme de formule voor de acceptatiekans α vereenvoudigt tot ¸ · π(tn ) ,1 . α = min π(sn−1 )

14


Empirische onderzoekingen leiden tot de aanbeveling de random-walk keuze zodanig te doen dat gemiddeld ongeveer 50% van de kandidaat-toestanden geweigerd worden. Voorbeeld 7.6. In dit voorbeeld lichten we de ‘continue’ versie van het MetropolisHastings algoritme toe aan de kansdichtheid π(x) die op een multiplicatieve constante na gegeven is door 1

1

x− 2 n e− 2 a/x

voor x > 0,

waarbij n = 5 en a = 4. Eerst beschrijven we hoe een stap van het algo2

ritme eruit ziet als we de onafhankelijkheidskeuze doen voor q(t) de χ1 dichtheid p 1 2/(πt)e− 2 t , t > 0 nemen. Bij de keuze s0 = 10 voor de begintoestand is de eerste stap: 2

(a) Een random trekking t1 uit de χ1 dichtheid wordt hgedaan, stel i dit geeft π(t1 )q(s0 ) t1 = 1.25. Dit resulteert voor de acceptatiekans α = min π(s0 )q(t1 ) , 1 in ´  1 2/(πs0 )e− 2 s0 ´ , 1 = 0.1987. α = min  ¡ ¢ ³p −2.5 −2/s0 − 12 t1 s0 e 2/(πt1 )e ¡

t−2.5 e−2/t1 1

¢ ³p

(b) Een random getal u tussen 0 en 1 wordt getrokken, stel u = 0.157. In dit geval is u ≤ α en wordt de kandidaat-toestand t1 geaccepteerd als nieuwe toestand s1 van de Markov keten. Het is leerzaam om het Metropolis-Hastings algoritme te proberen voor zowel de 2 keuze van de uniforme dichtheid op (0, 1000) als de keuze van de χ1 dichtheid voor de kandidaat-dichtheid. Deze uniforme dichtheid bedekt vrijwel de gehele kansmassa van de te simuleren dichtheid π( x), maar leidt tot een Markov keten met een heel slechte menging van de toestand doordat heel vaak de kandidaat2 toestand niet geaccepteerd wordt. De χ1 dichtheid geeft een betere menging maar alleen in een klein deel van het waardebereik van de dichtheid π(x): de 2 χ1 dichtheid heeft een veel te dunne staart in vergelijking met de dikke staart 2 van de dichtheid π(x) (de kans op een trekking groter dan 15 uit de χ1 dichtheid is 0.0001, terwijl de kansdichtheid π(x) nog een kansmassa 0.034 op (15, ∞) 2 heeft!). De tekortkomingen van de uniforme dichtheid en de χ1 dichtheid blijken ook uit schattingen voor de verwachtingswaarde van de dichtheid π(x) die na een miljoen trekkingen met het Metropolis-Hastings algoritme verkregen worden: deze schattingen zijn 3.75 en 3.04 terwijl de exacte waarde 4 is. Veel betere resultaten worden verkregen door voor de kandidaat-dichtheid q(x) een Paretodichtheid te kiezen met hetzelfde staartgedrag x−2.5 als de dichtheid π(x). De Pareto dichtheid wordt gegeven door (p/b)(b/x)p+1 voor x > b, waarbij de parameters b en p positieve getallen zijn. De passende keuze voor p is dus p = 1.5 en voor b is b = 1 een geschikte keuze. Als de stochastische variabele X hoort bij

15 Figuur 7.2 Menging in Metropolis-Hastings 100

100

90

90

80

80

70

70

60

60

50

50

40

40

30

30

20

20

10

10

0

0

100

200

300

400

500

0

0

100

200

300

400

500

de dichtheid π(x), dan wordt de kandidaat dichtheid q(t) = 1.5t−2.5 voor t > 1 gebruikt om de dichtheid πshif t (y) van de stochast Y = X +1 te simuleren, waarbij πshif t (y) = π(y − 1) voor y > 1. Deze aanpak geeft een Metropolis-Hastings algoritme met een goede menging van de toestanden in het gehele waardebereik van π(x) en leidt na een miljoen trekkingen tot de schatting 3.98 voor E(X). Het is uiterst instructief om voor elk van de gebruikte keuzes voor q(t) een plot te maken van zeg de eerste 500 waarnemingen uit de gesimuleerde Markov keten met s0 = 1 als begintoestand. In Figuur 7.2 geven we de op verticale as de waarden van de eerste 500 gesimuleerde toestanden zowel voor q(t) de uniforme dichtheid op (0, 1000) (linkerfiguur) als voor q(t) de Pareto dichtheid 1.5t−2.5 (rechterfiguur). Je ziet in de linkerfiguur dat de menging van de toestand heel slecht is voor de uniforme dichtheid: de kandidaat-toestand wordt vrijwel steeds verworpen. Voorbeeld 7.7. Stel je wilt een random element genereren uit een grote gecompliceerde ‘combinatorische’ verzameling V. Als het praktisch onuitvoerbaar is om dit rechtstreeks te doen, dan zou je dit kunnen doen met Metropolis-Hastings simulatie. Dit illustreren we voor de situatie dat V de verzameling is van alle Pn permutaties (s1 , . . . , sn ) van de getallen 1, . . . , n waarvoor j=1 jsj > a voor een gegeven constante a. Voor elk element s = (s1 , . . . , sn ) ∈ V, definiëren wij de omgevingsverzameling N (s) als alle elementen t = (t1 , . . . , tn ) ∈ V die uit s ontstaan door een verwisseling van twee posities. Bijvoorbeeld, bij n = 5 en a = 50 behoort (1, 2, 4, 3, 5) wel tot de omgevingsverzameling van (1, 2, 3, 4, 5) maar (5, 2, 3, 4, 1) niet. Pas nu het Metropolis-Hastings algoritme toe met q(t | s) =

1 |N (s)|

voor t ∈ N (s).

16


De evenwichtsverdeling die we door de toepassing van het algoritme willen bereiken is π(s) = 1/|V| voor s ∈ V. Dit betekent dat de acceptatiekans ¶ µ |N (s)| α(s, t) = min ,1 . |N (t)| Dus als de huidige toestand van de Markov keten s is, dan kies je random één van de buren uit N (s), stel het element t. Als t minder buren heeft dan s, dan blijf je in s. Anders, wordt een random getal U gegenereerd en is de volgende toestand van de Markov keten gelijk aan t als U < |N (s)|/|N (t)| en is gelijk aan s anders. Een random gekozen element uit V wordt verkregen na een voldoend groot aantal iteraties (de geconstrueerde Markov keten is aperiodiek).

1.0.3

De Gibbs sampler

De Gibbs sampler is een speciaal geval van het Metropolis-Hastings algoritme en wordt gebruikt om trekkingen te doen uit een multivariate kansdichtheid waarvan de univariate conditionele dichtheden bekend zijn. De Gibbs sampler vindt veelvuldig toepassing in de Bayesiaanse statistiek. Bij de Gibbs sampler ga je uit van een multivariate stochastische vector (X1 , . . . , Xd ) met simultane kansdichtheid π(x1 , . . . , xd ) = P (X1 = x1 , . . . , Xd = xd ). De univariate conditionele kansdichtheden πk (x|x1 , . . . , xk−1 , xk+1 , . . . , xd ) van de stochastische vector (X1 , . . . , Xd ) worden voor k = 1, . . . , d gedefinieerd door P (Xk = x | X1 = x1 , . . . , Xk−1 = xk−1 , Xk+1 = xk+1 , . . . , Xd = xd )

De aanname dat deze univariate conditionele dichtheden expliciet bepaald kunnen worden is essentieel voor de Gibbs sampler. Door uit de univariate conditionele dichtheden te trekken, genereert de Gibbs sampler een rij van opeenvolgende toestanden uit een Markov keten die π(x1 , . . . , xd ) als evenwichtsverdeling heeft. In de Gibbs sampler is de acceptatiekans altijd gelijk aan 1. Het algoritme is als volgt. Gibbs algoritme Stap 0. Kies een begintoestand x = (x1 , . . . , xd ). Stap 1. Genereer random een geheel getal k uit {1, . . . , d}. Doe een trekking y uit de univariate kansdichtheid πk (x|x1 , . . . , xk−1 , xk+1 , . . . , xd ). Laat y = (x1 , . . . , xk−1 , y, xk+1 , . . . , xd ). Stap 2. Laat x := y. Herhaal stap 1 met de nieuwe toestand x.

17 Dit is een toepassing van het Metropolis-Hastings algoritme met 1 q(y | x) = P (Xk = y | Xj = xj voor j = 1, . . . , d met j 6= k) d voor x = (x1 , . . . , xk−1 , xk , xk+1 , . . . , xd ) en y = (x1 , . . . , xk−1 , y, xk+1 , . . . , xd ). Voor de Gibbs sampler geldt dat de acceptatiekans α(x, y) altijd gelijk aan 1 is. Dit is eenvoudig in te zien. Op grond van q(y | x) = vinden we

1 π(y) d P (Xj = xj , j 6= k)

q(x | y) q(y | y)

=

π(x) π(y)

en

q(x | y) =

1 π(x) , d P (Xj = xj , j 6= k)

zodat ·

¸ π(y)π(x) α(x, y) = min , 1 = 1. π(x)π(y) Variant van het Gibbs algoritme Een veelgebruikte variant van bovenstaand algoritme is de Gibbs sampler waarin niet elke keer een random gekozen component van de toestandsvector aangepast wordt, maar in iedere iteratie achtereenvolgens elk van de componenten van de toestandsvector wordt aangepast. Dus als in de nde iteratie de toestandsvector (n) (n) (n) x(n) = (x1 , x2 , . . . , xd ) verkregen is, dan verloopt de (n + 1)ste iteratie als volgt: (n+1)

(n)

(n)

(n)

x1 is een random trekking uit π1 (x | x2 , x3 , . . . , xd ) (n+1) (n+1) (n) (n) x2 is een random trekking uit π2 (x | x1 , x 3 , . . . , xd ) . . . (n+1) (n+1) (n+1) (n+1) xd is een random trekking uit πd (x | x1 , x2 , . . . , xd−1 ). (n+1)

Dit geeft de aangepaste toestandsvector x(n+1) = (x1

(n+1)

, x2

(n+1)

, . . . , xd

) voor

de volgende iteratie.∗ Stel dat h(x1 , . . . , xd ) een gegeven functie is waarvoor we de verwachtingswaarde E[h(X1 , . . . , Xd )] willen schatten. Als de Gibbs de rij {x(k) , k = 1, 2, . . .} Pmsampler (k) van toestanden genereert, dan geeft (1/m) k=1 h(x ) een schatting voor de gezochte E[h(X1 , . . . , Xd )] wanneer m voldoende groot gekozen wordt. In het bijzonder stelt de Gibbs rij ons in staat om voor een individuele stochast Xk , zeg X1 , de marginale kansdichtheid of de verwachtingswaarde te schatten. Een ∗∗

Deze versie van de Gibbs Q sampler genereert een rij toestanden uit een Markov keten met d 1-stapsovergangskansen ps, t = j=1 πj (tj | t1 , . . . , tj−1 , sj+1 . . . , sd ).

18


na¨ıeve schatting voor de kansdichtheid π1 (x) van de stochast X1 wordt verkregen (k) door het kanshistogram te baseren op de waarden x1 voor k = 1, . . . , m uit de Gibbs rij {x(k) , k = 1, 2, . . . , m}. Een betere schatting wordt verkregen door gebruik te maken van de expliciete uitdrukking voor de univariate conditionele kansdichtheid π1 (x | x2 , x3 , . . . , xd ). De conditionele dichtheid stelt ons in staat π1 (x) te schatten door m

1 X (k) (k) πˆ1 (x) = π1 (x | x2 , . . . , xd ). m k=1 De theoretische achtergrond van deze schatter is π1 (x) = E[π1 (x | X2 , . . . , Xd )] (de wet van voorwaardelijke verwachtingen). De laatstgegeven schatter maakt gebruik van meer informatie en zal daardoor in het algemeen een betere schatter zijn met een lagere variantie. Voorbeeld 7.8. Veronderstel dat de simultane kansdichtheid π(x, y) van de stochasten X en Y op een multiplicative constante gegeven wordt door µ ¶ r x+α−1 y (1 − y)r−x+β−1 x

for x = 0, 1, . . . , r,

0 ≤ y ≤ 1,

waarbij r, α en β gegeven gehele positieve getallen zijn. Merk op dat de stochast X discreet verdeeld is en de stochast Y continu verdeeld. De Gibbs sampler werkt op dezelfde wijze voor discrete en continue verdelingen of mengsels van deze verdelingen. Stel nu dat we ge¨ınteresseerd zijn om de marginale kansdichtheid π1 (x) = P (X = x) van de stochast X te berekenen of de verwachting van X. Dit kunnen we doen met behulp van de Gibbs sampler. Het is vrij eenvoudig om uit de vorm van de simultane kansdichtheid π(x, y) af te leiden dat π1 (x | y) is de binomiale kansdichtheid met parameters r en y π2 (y | x) is de beta dichtheid met parameters x + α en r − x + β. Een lange Gibbs rij (x0 , y0 ), (x1 , y1 ), . . . , (xm , ym ) voor een voldoende grote m wordt als volgt gegenereerd. Kies een geheeltallige beginwaarde x0 tussen 0 en r. De rest van de rij wordt iteratief verkregen door afwisselend een waarde yj te trekken uit de beta dichtheid π2 (y | xj ) en een waarde xj+1 uit de binomiale dichtheid π1 (x | yj ). Codes om uit de binomiale dichtheid en de beta dichtheid random trekkingen te doen zijn alom beschikbaar. Voor de getalwaarden r = 16, α = 2, en β = 4, hebben we de tweede versie van de Gibbs sampler gebruikt om de marginale dichtheid π1 (x) te schatten. In totaal zijn m = 250000 waarnemingen voor de toestandsvector (x, y) gegenereerd. In Figuur 7.3 geven we het gesimuleerde kanshistogram voor de marginale dichtheid π1 (x), waarbij we

1.1. OPGAVEN

19

Figuur 7.3 De gesimuleerde en exacte kanshistogrammen 0.100 0.075 0.050 0.025

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

0.100 0.075 0.050 0.025

de conditionele schatter π ˆ1 (x) ¡gebruikt hebben. In dit voorbeeld wordt πˆ1 (x) ¢ x P m r 1 gegeven door πˆ1 (x) = m k=1 x yk (1 − yk )r−x voor x = 0, 1 . . . , r. Onder het gesimuleerde kanshistogram geven we ter vergelijking het histogram met de exacte waarden van de π1 (x) (in het beschouwde voorbeeld is het mogelijk de exacte ¡¢ voor x = 0, 1, . . . , r). kansen te berekenen uit π1 (x) = xr (α+β−1)!(x+α−1)!(r−x+β−1)! (α−1)!(β−1)!(α+β+r−1)!

1.1

Opgaven

1. Simuleer met het Metropolis-Hastings algoritme het kanshistogram van π(1) = 0.2 en π(2) = 0.8 door het gooien met een zuivere munt (d.w.z. q(s|t) = 0.5 voor s, t=1, 2) en zie hoe snel convergentie plaatsvindt. 2. Veronderstel dat de simultane kansdichtheid van de continue stochasten X1 en X2 op een multiplicative constante na gegeven wordt door 1

2 2

2

2

e− 2 (x1 x2 +x1 +x2 −7x1 −7x2 )

voor − ∞ < x1 , x2 < ∞.

Experimenteer met Markov keten Monte Carlo simulatie om de verwachtingswaarde, de spreiding en de marginale kansdichtheid van X1 te vinden. Pas zowel de Gibbs sampler toe als het Metropolis-Hastings algoritme met de random-walk keuze (x1 , x2 ) + (Z1 , Z2 ) met Z1 en Z2 onafhankelijke N (0, a2 ) verdeelde stochasten. Probeer verschillende waarden van a (zeg, a = 0.02, 0.2, 1 en 5) en zie hoe de menging van de toestand verloopt en wat de gemiddelde waarde van de acceptatiekans is. 2. In een actuarieel model hebben de stochasten X, Y en N een simultane kansdichtheid π(x, y, n) die op een multiplicatieve constante na gegeven wordt door

20

CHAPTER 1. MARKOV KETEN MONTE CARLO METHODEN µ ¶ n x+α−1 λn y (1 − y)n−x+β−1 e−λ x n!

voor x = 0, 1, . . . , n, 0 < y < 1 en n = 0, 1, . . .. De stochast N representeert het aantallen polissen in een portefeuille, de stochast Y representeert de kans dat een gegeven polis tot een claim leidt (elke polis heeft dezelfde claimkans en de polissen gedragen zich onafhankelijk van elkaar), en de stochast X representeert het aantal claims dat resulteert van de polissen in de portefeuille. Voor de data α = 2, β = 8 en λ = 50 simuleer met de Gibbs sampler het kanshistogram van de marginale dichtheid van X alsmede de verwachting en spreiding van X zowel met een na¨ıeve schattingsprocedure als met de schattingsprocedure die de expliciete formule voor de univariate conditionele dichtheid π1 (x | y, n) gebruikt. Ga daartoe eerst na dat de univariate conditionele kansdichtheden van X, Y en N worden gegeven door: π1 (x | y, n) is een binomiale(n, y) dichtheid, π2 (y | x, n) is een beta(x + α, n − x + β) dichtheid en π3 (n | x, y) is een Poisson(λ(1 − y)) dichtheid verschoven naar x.

Chapter 2 Stochastic Dynamic Programming 2.1

Stochastische dynamische programmering

De toepassing van dynamische programmering op stochastische sequentiële beslissingsproblemen is conceptueel niet moeilijker dan in het geval van deterministische sequentiële beslissingsproblemen. Bij problemen met deterministische toestandsovergangen kan je van tevoren al bepalen welke reeks van beslissingen je zult nemen. Bij stochastische dynamische programmeringsproblemen zijn de toekomstige toestanden onzeker en is het dus niet mogelijk om reeds van tevoren te zeggen welke reeks van beslissingen genomen zal worden. Wat we nodig hebben is een verzameling van conditionele beslissingen van de vorm ‘Als de gerealiseerde toestand gelijk is aan . . . , neem dan de beslissing . . . .’ Een dergelijke strategie kan verkregen worden door de achterwaartse recursie van dynamische programmering toe te passen. De redenering achter de recursieve relatie voor stochastische dynamische programmeringsproblemen is in essentie dezelfde als voor deterministische problemen. De beste manier om dit te laten zien is met behulp van enkele karakteristieke voorbeelden van stochastische dynamische programmering.

2.1.1

Een dobbelspel en optimaal stoppen

Een interessant spel is het volgende. Het spel bestaat uit het maximaal zes keer gooien van een zuivere dobbelsteen. Na elke worp moet je beslissen of je doorgaat of stopt. De uitbetaling van het spel is het aantal geworpen ogen van de laatste worp voordat je stopte. Welke strategie moet je volgen om de verwachte uitbetaling te maximaliseren? Het probleem is een typisch voorbeeld van een stochastisch sequentieel beslissingsprobleem. De toestand van een dynamisch proces wordt geobserveerd op discrete tijdstippen. Na elke observatie van de toestand wordt een beslissing genomen. Dan wordt een directe opbrengst ontvangen, die alleen gebaseerd 21

22

CHAPTER 2. STOCHASTIC DYNAMIC PROGRAMMING

is op de toestand op dat moment en de genomen beslissing. Vervolgens gaat het proces naar de volgende toestand overeenkomstig een gegeven kansverdeling. Welke strategie moet gekozen worden opdat de verwachte waarde van de totaal te verkrijgen opbrengst over een gegeven periode met eindige lengte maximaal is? Dynamische programmering stelt je in staat om dergelijke problemen op te lossen met gebruik van een recursief algoritme. Op dezelfde manier als in de vorige paragrafen voor deterministische sequentiële beslissingsproblemen, wordt het oorspronkelijke probleem opgedeeld in een reeks van geneste subproblemen die recursief aan elkaar gekoppeld worden. Om dat te doen hebben we het basisbegrip waardefunctie nodig. Voor het beschouwde probleem definiëren wij voor k = 0, 1, . . . , 6 de functie fk (i) door fk (i) = het maximale verwachte aantal punten als nog k worpen te gaan zijn en het aantal punten van de laatste worp gelijk aan i is. Het doel is om f6 (0) en de optimale strategie te bepalen. Om de recursieve relatie voor fk (i) te bepalen, redeneer je als volgt. Er zijn twee mogelijke acties ‘stop’ en ‘ga door’ in de huidige beslissingssituatie waarin nog k worpen te gaan zijn en i punten verkregen zijn in de vorige worp. Als de beslissing om te stoppen genomen wordt, dan wordt een onmiddellijke beloning van i ontvangen en is het spel afgelopen. Als wij doorgaan met het spel, dan wordt er nog geen beloning ontvangen en zal de volgende toestand van het proces gelijk aan j zijn met kans 1/6 voor j = 1, . . . , 6. Als voor doorgaan wordt gekozen en vervolgens optimaal verder wordt gespeeld in deP resterende k−1 worpen, dan is de verwachte uitbetaling voor het spel gelijk aan 6j=1 fk−1 (j)/6. Dus vinden we de volgende recursieve relatie: ½ ¾ 6 1X fk (i) = max i, fk−1 (j) . 6 j=1 Startend met f0 (i) = i voor alle i, wordt achtereenvolgens voor k = 1, . . . , 5 de waardefunctie fk (i) berekend voor i = 1, . . . , 6. Tenslotte volgt de verwachte uitbetaling van een spel uit 6

f6 (0) =

1X f5 (j). 6 j=1

Voor elke combinatie (k, i) is het nodig om de optimale beslissing, zeg dk (i), waarvoor het maximum in de recursieve relatie voor fk (i) wordt bereikt te bewaren. Dit ‘handboek’ van beslissingen geeft je een optimale strategie. Numerieke berekeningen Het is leerzaam om de berekeningen voor dit specifieke probleem uit te voeren. De algoritme wordt gestart met f0 (i) = i voor i = 1, . . . , 6, waarbij d0 (i) =‘stop’ voor

2.1. STOCHASTISCHE DYNAMISCHE PROGRAMMERING

23

alle i. De waardefunctie fk (i) wordt achtereenvolgens berekend voor k = 1, . . . , 5: ½ f1 (i) = ½ f2 (i) = ½ f3 (i) = ½ f4 (i) = ½ f5 (i) =

3.5 i

voor i = 1, 2, 3 voor i = 4, 5, 6

met d1 (i) = ‘ga door’, met d1 (i) = ‘stop’,

4.25 i

voor i = 1, 2, 3, 4 voor i = 5, 6


4.667 i

voor i = 1, 2, 3, 4 voor i = 5, 6


4.944 i

voor i = 1, 2, 3, 4 voor i = 5, 6


5.130 i

voor i = 1, 2, 3, 4, 5 met d5 (i) = ‘ga door’, voor i = 6 met d5 (i) = ‘stop’.

Tenslotte bereken je de maximale verwachte uitbetaling als f6 (0) = 5.275. Dit heeft tot gevolg dat het gunstig is om dit spel te spelen als het spel een inzet s heeft met 0 < s < 5.275. Dan zal, op grond van de wet van de grote aantallen, je feitelijke gemiddelde winst per spel willekeurig dicht bij 5.275 − s komen als het spel voldoende vaak wordt gespeeld en de optimale strategie wordt gebruikt. De berekeningen van het dynamisch-programmeringsalgoritme laten zien dat de optimale strategie een eenvoudige structuur heeft en gekarakteriseerd wordt door de kengetallen s1 = 4, s2 = 5, s3 = 5, s4 = 5, en s5 = 6. Als nog m worpen te gaan zijn, dan schrijft de strategie voor te stoppen als de vorige worp sm of meer punten heeft opgeleverd en door te gaan anders (m = 1, . . . , 5). Ter afsluiting van deze subparagraaf beschouwen we nog enkele varianten van het behandelde spel. In de eerste variant bestaat het spel uit het ten hoogste M keer werpen van een onzuivere dobbelsteen, waarbij M een gegeven geheel getal is. Elke worp met de dobbelsteen levert j punten op met een gegeven kans p(j) voor j = 1, . . . , 6. Als we de waardefunctie definiëren als hierboven, dan volgt de recursieve relatie ½ fk (i) = max i,

6 X

¾ fk−1 (j) p(j)

j=1

voor k = 1, . . . , M met de randconditie f0 (i) = i voor alle i. Het wordt nu aan jezelf over gelaten om de dynamische-programmeringsformulering te vinden voor de tweede variant, waarbij twee zuivere dobbelstenen worden gegooid en de som van de twee dobbelstenen telt.

24


2.1.2

Het spel rood en zwart

Stel dat je naar een casino gaat omdat je een bepaalde som geld nodig hebt voor de volgende ochtend. Gokken is de laatste mogelijkheid om het geld te verkrijgen. Je besluit om het spel met rood en zwart te gaan spelen. Vanwege de tijdslimiet kan je nog maar een eindig aantal keren, zeg n keer, spelen. Elke keer kan je elk geheel bedrag in euro’s inzetten tot aan je vermogen. Je wint je inzet plus je eigen geld terug met een gegeven kans p en je verliest je inzet met kans q = 1 − p. Je oorspronkelijke kapitaal is gelijk aan A euro en je doel is om tenminste B euro te bereiken met B > A. Wat is de optimale inzetstrategie als je de kans wilt maximaliseren op het verkrijgen van het benodigde geld in niet meer dan n weddenschappen? Dit probleem kan ge¨ınterpreteerd worden als een stochastisch sequentieel beslissingsprobleem. De toestand is gelijk aan i als je huidige vermogen gelijk is aan i euro. Voor toestand i zijn de mogelijke beslissingen d = 0, 1, . . . , i, waarbij beslissing d correspondeert met het inzetten van d euro. Om dit probleem met behulp van dynamische programmering op te lossen, definiëren we voor k = 1, . . . , n de waardefunctie fk (i) = de maximale kans op het bereiken van een vermogen van tenminste B euro als je nog k keer mag inzetten en je huidige vermogen gelijk is aan i euro. De volgende redenering wordt gevolgd om een recursieve relatie voor fk (i) te vinden. Stel dat je d euro inzet met k weddenschappen te gaan en vervolgens optimaal inzet in de resterende k − 1 weddenschappen. Onder de conditie dat de volgende toestand gelijk is aan s, bereik je je doel met kans fk−1 (s). De volgende toestand is of i + d of i − d met respectievelijke kansen p en q. Dus in de situatie dat er k weddenschappen te gaan zijn en je huidige vermogen gelijk is aan i euro, bereik je je uiteindelijke doel met kans pfk−1 (i + d) + qfk−1 (i − d) als je d euro inzet en vervolgens optimaal verder speelt in de resterende k − 1 spelen. Het maximaliseren van deze kans over alle mogelijke d geeft de recursie fk (i) = max {pfk−1 (i + d) + qfk−1 (i − d)}. d=0,1,...,i

Elke beslissing d waarvoor het maximum in deze recursieve relatie wordt bereikt is een optimale beslissing voor de situatie waarin er nog k weddenschappen te gaan zijn en je huidige vermogen gelijk is aan i euro. Door op recursieve wijze de waardefunctie te berekenen, verkrijg je een optimale inzetstrategie. De recursie wordt gestart met f0 (i) = 1 voor i ≥ B en f (i) = 0 voor i < B. Alles-of-niets strategie Als de winstkans p voldoet aan p ≤ 12 , dan geldt dat de optimale strategie is om domweg je hele vermogen i in te zetten als i < B/2 en B − i als je huidig

2.2. INVESTERINGSPROBLEEM EN DE KELLY STRATEGIE

25

vermogen i groter is dan of gelijk aan B/2. Het bewijs van dit resultaat vereist zeer diepgaande wiskunde. Een intu¨ıtieve verklaring voor de optimaliteit van de alles-of-niets strategie voor het geval van een ongunstig spel is dat deze je geld zo kort mogelijk aan het huisvoordeel van het casino blootstelt en op die manier je winstkans zo hoog mogelijk maakt. Numeriek voorbeeld Beschouw ter illustratie het getallenvoorbeeld A = 1, B = 5, p = 0.5, n = 7. Tabel 5.6 Numerieke resultaten k/i 1 2 3 4 5 6 7

1 0 0 0.125 0.1875 0.1875 0.1875 0.1953

fk (i) 2 3 0 0.5 0.25 0.5 0.375 0.5 0.375 0.5625 0.375 0.5938 0.3906 0.5938 0.3984 0.5938

4 0.5 0.75 0.75 0.75 0.7813 0.7969 0.7969

1 0-1 0-1 1 1 0-1 0-1 1

dk (i) 2 3 4 0-1-2 2-3 1-2-3-4 1-2 2-3 1 2 1-2-3 0-1 0-2 1-2 0-1 1-2 2 1 1-2 2 1 2 0-1-2 1

In tabel 5.6 zijn de berekeningen samengevat. De tabel geeft zowel de waardefunctie fk (i) als de optimale beslissing dk (i). Zoals uit de resultaten blijkt zijn er verscheidene optimale inzetstrategieën, waaronder de alles-of-niets strategie. De maximale kans om de benodigde 5 euro te verkrijgen is 0.1953 als je start met 1 euro en er niet meer dan 7 weddenschappen zijn toegestaan.

2.2

Investeringsprobleem en de Kelly strategie

Als inleiding tot het investeringsprobleem, beschouwen we het volgende scenario. Tijdens de internethausse op de beurs aan het eind van de vorige eeuw gingen talloze dotcombedrijven naar de beurs. Vooral in de eerste week kon de beurskoers van een nieuw dotcombedrijf sterk fluctueren. Stel eens dat gemiddeld van de helft van de nieuwe dotcombedrijven op de beurs de koers in de eerste week met 80% stijgt en van de andere helft met 60% daalt. Dit betekent dat elke dollar ge¨ınvesteerd in een nieuw dotcombedrag na één week 0.5 × $1.8 1 + 0.5 × $0.4 = $1.1 als verwachte waarde heeft, een verwachte winstsstijging van 10%. Je hebt een beginkapitaal van $10 000 om te investeren in dit soort bedrijven. Stel eens dat je de volgende strategie gaat aanhouden voor de komende 52 weken: aan het begin van elke week beleg je het gehele huidige kapitaal in een nieuw dotcombedrijf en aan het eind van de week verkoop je het aandeel weer. Wat is de meest waarschijnlijke waarde van je kapitaal na 52 weken? De lezer wordt uitgenodigd om op dit moment niet verder te lezen en voor zichzelf een schatting te maken. Velen denken dat de meest waarschijnlijke waarde in de buurt van 15 duizend dollar zal liggen en achten de kans vrijwel gelijk aan nul dat je na

26


52 weken minder dan je beginkapitaal van 10 duizend dollar hebt. De meest waarschijnlijke waarde van je kapitaal na 52 weten is echter $1.95! Dit kun je direct inzien door te bedenken dat het aantal stijgingen in 52 weken binomiaal verdeeld is met parameters n = 52 en p = 0.5. Dus de meest waarschijnlijke waarde van je kapitaal na 52 weken is (1.8)26 (0.4)26 × $10000 = $1.95. Als je elke keer je gehele kaptaal herinvesteer, dan is de kans op een eindkapitaal van niet meer dan $1.95 gelijk aan de binomiaalkans 0.555 en is de kans op een kapitaal hoger dan je beginkapitaal van 10 duizend dollar gelijk aan 0.058. De winststijging van gemiddeld 10% per week is misleidend: de factor 1.8 × 0.4 is bepalend en deze factor is kleiner dan 1. Bij investeren over een langere periode is het echter optimaal om elke keer eenzelfde vaste fractie van je kapitaal in te zetten. In het voorgaande voorbeeld is het optimaal om elke keer dezelfde 5 fractie α∗ = 24 van je kapitaal in te zetten. In dat geval kan worden berekend dat de kans op een kapitaal van niet meer dan $1.95 nihil is, teerwijl de kans op 5 een eindkapitaal hoger dan 10 duizend dollar bijna 70% is. De waarde α∗ = 24 volgt uit een algemene formule die bekend staat onder de naam Kelly formule. Deze formule heeft betrekking op de situatie dat je een reeks van onafhankelijke investeringsmogelijkheden hebt, waarbij je bij elke investering het ge¨ınvesteerde bedrag f1 keer terugkrijgt met kans p en f2 keer met kans 1 − p met waarden voor p, f1 en f2 die voldoen aan f1 > 1,

0 ≤ f2 < 1

en

pf1 + (1 − p)f2 > 1.

In het bovenstaande voorbeeld is p = 0.5, f1 = 1.8 en f2 = 0.4. Bij toepassen van de Kelly strategie is de groei van je kapitaal op de lange duur maximaal en overschrijdt je kapitaal uiteindelijk elke waarde. Dit resultaat kunnen we plausibel maken door de situatie te beschouwen dat we een reeks van N investeringsmogelijkheden hebben en het doel is om de verwachtingswaarde van de logaritme van het eindkapitaal te maximaliseren. De keuze voor een logaritmische utiliteitsfunctie wordt niet verder toegelicht, maar is te verklaren op grond van economische overwegingen. Definieer de waardefunctie Lk (x) = de maximale verwachte utiliteitswaarde van je eindkapitaal als er nog k keer ge¨ınvesteerd mag worden en je huidige kapitaal gelijk is aan x. voor k = 0, . . . , N en x > 0, waarbij L0 (x) = ln(x). Wij veronderstellen dat het kapitaal ‘oneindig deelbaar’ is, d.w.z. elke fractie van je huidige kapitaal mag elke keer ge¨ınvesteerd worden. De beslissing in elke toestand representeren we als de fractie van het huidige kapitaal dat je investeert. We hebben nu te doen met een dynamisch-programmeringsprobleem waarvoor de verzameling van toestanden een ‘continue’ verzameling is evenals de verzameling van mogelijke beslissingen. Voor het principe van dynamische programmering maakt dit echter niets uit. Het volgende resultaat zullen we nu bewijzen. Om een recursieve relatie


27

te vinden voor Lk (x) redeneren wij als gebruikelijk. Stel dat je een fractie α inzet van je huidige kapitaal x en vervolgens optimaal handelt in de resterende k − 1 investeringsmogelijkheden. Je vermogen als nog k − 1 perioden te gaan zijn, is met kans p gelijk aan x−αx+f1 ×αx en met kans 1−p gelijk aan x−αx+f2 ×αx zodat de verwachte waarde van de logaritme van je eindvermogen gegeven wordt door pLk−1 ([1 + (f1 − 1)α]x) + (1 − p)Lk−1 ([1 + (f2 − 1)α]x) . Maximaliseren we deze uitdrukking over α, dan vinden we de recursie Lk (x) = max {pLk−1 ([1 + (f1 − 1)α]x) + (1 − p)Lk−1 ([1 + (f2 − 1)α]x)}. 0≤α≤1

Het volgende resultaat zullen we nu bewijzen. Voor alle k = 1, . . . , N en x > 0 geldt dat de maximaliserende beslissing α∗ in de optimaliteitsvergelijking en de waardefunctie Lk (x) gegeven worden door µ ∗

α = min

¶ pf1 + (1 − p)f2 − 1 ,1 (f1 − 1)(1 − f2 )

en

Lk (x) = kc + ln(x),

waarbij c = p ln(1 + (f1 − 1)α∗ ) + (1 − p) ln(1 + (f2 − 1)α∗ ) > 0. Met andere woorden, de optimale strategie is altijd dezelfde fractie a∗ van je huidige kapitaal in te zetten, ongeacht de hoogte van je kapitaal en het aantal investeringsmogelijkheden dat nog te gaan is. Een waarlijk opmerkelijk resultaat! Het bewijs van dit resultaat gaat met behulp van inductie. Het inductiebewijs start met n = 1. Aangezien L0 (x) = ln(x) en ln(ab) = ln(a) + ln(b), kunnen wij de uitdrukking voor L1 (x) schrijven als L1 (x) = ln(x) + max {p ln(1 + (f1 − 1)α) + (1 − p) ln(1 + (f2 − 1)α)}. 0≤α≤1

Beschouw de functie h(α) = p ln(1 + (f1 − 1)α) + (1 − p) ln(1 + (f2 − 1)α) voor 2 −1 0 ≤ α ≤ 1. De oplossing van h0 (α) = 0 is α0 = pf(f11+(1−p)f . De function h(α) −1)(1−f2 ) is concaaf in α met h(0) = 0. Dit betekent dat de functie h(α) op (0, 1) het maximum aanneemt in α∗ = min(α0 , 1). Verder zien we dat L1 (x) = ln(x) + c. Veronderstel dat de bewering bewezen voor n = 1, . . . , k − 1. Substituteren we Lk−1 (x) = (k − 1)c + ln(x) in de recursie vergelijking voor Lk (x), dan volgt dat

Lk (x) = ln(x) + (k − 1)c + max {p ln(1 + (f1 − 1)α) + (1 − p) ln(1 + (f2 − 1)α)}. 0≤α≤1

De vergelijking voor Lk (x) is dezelfde als die voor L1 (x), op de additieve constante (k − 1)c na. Dit heeft tot gevolg dat de maximaliserende waarde α in de recursieve vergelijking voor Lk (x) ook gegeven wordt door α = α∗ en dat Lk (x) = (k − 1)c + ln(x) + c = kc + ln(x). Hiermee is het bewijs rond.

28


2.2.1

Dynamische programmering met een stochastische planningsduur

Een populair dobbelspel in Amerika is de zogenoemde ‘Game of Pig’. Dit spel wordt gespeeld door twee personen die om de beurt aan zet zijn. De speler die aan zet is, gooit elke keer één dobbelsteen tegelijk maar kan de dobbelsteen meerdere malen werpen in zijn speelbeurt. Als de dobbelsteen de uitkomst 1 geeft, dan is de speelbeurt van de speler over en is de andere speler aan de beurt, anders moet de speler beslissen de dobbelsteen nogmaals te gooien of te stoppen. Ingeval de speler een speelbeurt moet stoppen door het gooien van een 1, dan vervallen de punten die in de speelbeurt gegooid zijn en blijft zijn totaalscore hetzelfde als aan het begin van de speelbeurt. Ingeval de speler besluit de speelbeurt te stoppen terwijl geen 1 gegooid is, dan wordt het totale aantal punten gegooid in de speelbeurt aan de totaalscore van de speler toegevoegd. In het twee-spelers geval is de speler die voor eerst 100 punten bereikt de winnaar. In dit voorbeeld beschouwen we de 1-speler versie en zijn we ge¨ınteresseerd in het minimale verwachte aantal speelbeurten nodig om 100 punten te bereiken. Dit probleem is een interessant en leerzaam voorbeeld van een dynamisch programmeringsprobleem met een planningsduur die stochastisch is en die eindigt bij het bereiken van een absorberende toestand. Dit probleem is een interessant en leerzaam voorbeeld van een dynamisch programmeringsprobleem met een planningsduur die stochastisch is en die eindigt bij het bereiken van een absorberende toestand. We zeggen dat de toestand (i, 0) is als je voor het begin van een nieuwe speelbeurt staat en je nog i punten nodig hebt om je doel van G = 100 punten te bereiken (je huidige totaalscore is dus 100−i punten). Voor de situatie dat een speelbeurt gaande is, wordt de toestand van het systeem wordt gedefinieerd door s = (i, k) met i het aantal punten dat je aan het begin van de huidige speelbeurt nog nodig had voor je einddoel en k (≥ 2) het aantal punten tot nu toe gegooid in de huidige beurt zonder een 1 gegooid hebben in die beurt. In elke toestand (i, k) zijn er twee beslissingen ‘stoppen’ of ‘gooien’, waarbij uiteraard de beslissing ‘gooien’ de enige mogelijke beslissing is in toestand (i, 0). De waardefunctie V (s) wordt gedefinieerd door V (s) = de minimale verwachte aantal nog nieuw te starten speelbeurten om 100 punten te bereiken vanuit toestand s. Het doel is V (100, 0) te vinden met de bijbehorende optimale strategie. De optimaliteitsvergelijking voor V (s) wordt verkregen met het gebruikelijke argument in dynamische programmering. Voor toestand (i, k) met k = 0. 6

X1 1 V (i, r). V (i, 0) = 1 + V (i, 0) + 6 6 r=2


29

Voor k ≥ 2 en i − k > 0, 6

V (i, k) = min[V (i − k, 0),

X1 1 V (i, 0) + V (i, k + r)], 6 6 r=2

waarbij V (i, l) = 0

voor (i, l) met i − l ≤ 0.

De eerste term in het rechterlid van de optimaliteitsvergelijking voor V (i, k) behoort bij de beslissing ‘stoppen’ in toestand (i, k) en daarna optimaal verder handelen vanuit toestand (i − k, 0) en de tweede term behoort bij de beslissing ‘gooien’ in toestand (i, k) en daarna optimaal verder handelen vanuit de resulterende toestand. Een nadere beschouwing van de optimaliteitsvergelijking voor V (s) leert dat deze vergelijking geen recursieve vergelijking is door het ontbreken van een tijdsindex die elke keer met 1 verandert. In de dynamische programmering is er een standaardmethode om een optimaliteitsvergelijking als bovenstaande op te lossen. Dit is de methode van successieve substitutie. Je begint met een willekeurig gekozen functie V0 (s), meestal kies je V0 (s) = 0 voor alle s. Vervolgens worden de functies V1 (s), V2 (s), . . . recursief berekend via 6

X1 1 Vn (i, 0) = 1 + Vn−1 (i, 0) + Vn−1 (i, r), 6 6 r=2 en 6

X1 1 Vn (i, k) = min[Vn−1 (i − k, 0), V (i, 0) + Vn−1 (i, k + r)]. 6 6 r=2 Op grond van een basisresultaat uit de theorie van de dynamische programmering geldt lim Vn (s) = V (s) voor alle s. n→∞ P In de praktijk stop je de berekeningen zodra s |Vm (s) − Vm−1 (s)| klein genoeg is, zeg kleiner dan 0.001 (merk op dat je in feite te maken hebt met een genest stelsel optimaliteitsvergelijkingen die achtereenvolgens voor i = 1, 2, . . . , 100 elk apart kunnen worden opgelost). Voor de minimale verwachte waarde van het aantal speelbeurten om 100 punten te bereiken, vinden we V (100, 0) = 12.545. De bijbehorende optimale strategie is tamelijk complex. Verrassend is dat in dit probleem een simpele heuristiek te geven die bijna even goed is als de optimale strategie. De heuristiek is om een speelbeurt te stoppen zodra 20 of meer punten in de beurt behaald zijn. Voor de ‘stop-bij-20’ heuristiek is het verwachte aantal speelbeurten gelijk aan 12.637, een relatief verschil van slechts 0.7% met de minimale waarde.∗ De ‘stop-bij-20’ regel leid je af met een klassiek recept voor heuristieken: doe net alsof nog één stap te gaan is. Stel dat je in de huidige ∗

Definieer je voor de ‘stop-bij-20’ heuristiek H(s) als de verwachte waarde van het

30


beurt k punten hebt en nogmaals gooit, dan is de directe verwachte toename in punten gelijk aan 56 × 4 en is het verwachte aantal punten dat je verspeelt gelijk aan 16 × k; de eerste k waarvoor 61 × k ≥ 56 × 4 is k = 20.

2.3

Opgaven

1. In de populaire TV-quiz ‘Eerst denken, dan doen’ word je verzocht een aantal vragen te beantwoorden. Er zijn ten hoogste K vragen. Op basis van je ervaring schat je dat de kans op het correct beantwoorden van vraag k gelijk is aan pk voor k = 1, . . . , K. De quiz hanteert het principe ‘alles of niets’. Een fout antwoord betekent dat je niet langer meedoet en niets wint. Als je vraag k correct beantwoordt, heb je de keus tussen stoppen en ak duizend euro ontvangen of doorgaan met de volgende vraag. Welke strategie moet je volgen om je verwachte opbrengst te maximaliseren? Formuleer een dynamisch-programmeringsalgoritme. Los dit probleem op voor de getalwaarden pk = 0.95, 0.90, 0.85, 0.75, 0.50, 0.70 en 0.75 en ak = 5, 15, 25, 35, 50, 75 en 100 voor k = 1, . . . , 7. 2. Veronderstel in opgave 1 dat je nu één vraag mag missen. Als je stopt na vraag k, dan krijg je ak of ak−1 duizend euro afhankelijk van of je tot dan geen enkele of één vraag gemist hebt. Los het probleem opnieuw op. 3. Stel dat je een aantal investeringsmogelijkheden aangeboden krijgt. Bij elke gelegenheid kan je investeren in ten hoogste een van de twee projecten A en B. Een investering in het risicovolle project A geeft een winst van 100% met kans pA en gaat verloren met kans 1 − pA . Je kan nooit geld verliezen als je in project B investeert. Een investering in project B levert een winst op van 100% met kans pB (< pA ) en de investering krijg je terug zonder winst met kans 1 − pB . Er zijn niet meer dan N investeringsmogelijkheden. Bij elke gelegenheid kan je elk geheeltallig bedrag in euro’s investeren tot en met je vermogen op dat moment, waarbij je winsten van eerdere investeringen mag herinvesteren. Welke investeringsstrategie moet je volgen opdat de kans maximaal is op het verkrijgen van S euro, als je huidige vermogen gelijk is aan R euro ? Formuleer een dynamisch-programmeringsalgoritme. Los dit probleem op voor de numerieke gegevens pA = 0.5, pB = 0.3, N = 5, R = 1 en S = 5. 4. (a) Stel dat je een randomgenerator gebruikt om elke keer aselect een getal uit de gehele getallen 1, . . . , 1000 te trekken. Na elke trekking heb je de optie om te stoppen of door te gaan met de volgende trekking. Je winst is gelijk aan het getal dat je in de laatste trekking hebt getrokken voordat je stopt. Gebruik dynamische programmering om de maximale verwachte waarde van je winst te resterende aantal benodigde speelbeurten inclusief de huidige speelbeurt om 100 punten te bereiken vanuit toestand s, dan kan de functie H(s) op vergelijkbare wijze berekend worden als de waardefunctie V (s) (in plaats van het minimum van twee termen heb je in de functionaalvergelijking voor H((i, k)) slechts één term omdat in elke toestand (i, k) slechts één beslissing mogelijk is). Dit leidt tot H(100, 0) = 12.637.

2.3. OPGAVEN

31

bepalen als je maximaal 20 trekkingen mag doen. (b) Een spel bestaat uit het ten hoogste zes maal werpen van twee zuivere dobbelstenen. Wat is een optimale stopregel als de uitbetaling gelijk is aan het aantal ogen dat geworpen is in de laatste worp? 5. Stel dat je een optie bezit om één aandeel van een bepaald fonds te kopen tegen een vaste prijs c gedurende een van de komende N dagen. Je bent niet verplicht om de optie te gebruiken. Je maakt een winst van s−c als je het aandeel koopt op een dag dat de koers gelijk is aan s. De huidige prijs van het aandeel is s0 en de koersen veranderen elke dag met 1, 0 of −1 met respectievelijke kansen p, q en r met p + q + r = 1. Welke strategie maximaliseert je verwachte winst? Formuleer een dynamisch-programmeringsalgoritme. Los het probleem op voor de numerieke gegevens c = 10, N = 7, s0 = 10, p = 0.3, q = 0.2 en r = 0.5. 6. Je mag tot K keer een zuivere munt opgooien, waarbij K vooraf gegeven is. Na elke worp moet je besluiten te stoppen of door te gaan. Je uitbetaling wordt in de eenheid van 1000 euro gegeven door de fractie worpen met kop op het moment dat je stopt. Bereken voor zowel K = 10, 25, 50 als 100 de maximale verwachte opbrengst en een optimale strategie. 7. Bij de verkoop van je huis wordt je geconfronteerd met de situatie dat je achtereenvolgens biedingen 1, 2, . . . , N op je huis zal krijgen met N een vastgelegde waarde. Als een bieding gedaan is, moet je direct beslissen of je deze accepteert of niet (bij de allerlaatste bieding heb je geen keus, die moet altijd geaccepteerd worden). Je kunt niet later terugkomen op een bod dat je eerder afgewezen hebt. De hoogtes van de biedingen zijn onafhankelijk van elkaar en elke bieding is een trekking uit een gegeven kansverdeling {aj }. • Formuleer een DP-algoritme waarmee berekend kan worden hoe je te werk moet gaan als je de verwachtingswaarde van het geaccepteerde bod wil maximaliseren de optimale strategie kan worden berekend. • Veronderstel nu dat je doel is om maximale kans het allerhoogste bod dat in de N biedingen zou optreden in de wacht te slepen. Wat wordt nu de waardefunctie en het DP-algoritme? 8. Beschouw opnieuw de 1-speler versie van de ‘Game of Pig’ uit paragraaf 5.8.4. Bereken met dynamische programmering voor N = 7, 10, 15 en 20 de maximale kans om 100 punten te bereiken als de speler over niet meer dan N speelbeurten de beschikking heeft. 9. Beschouw de volgende variant van de ‘Game of Pig’ uit paragraaf 5.8.4. Een speler gooit elke keer met twee dobbelstenen tegelijk, waarbij in een speelbeurt ook weer meerdere malen geworpen kan worden. Als een speler een 1 gooit, dan eindigt de beurt met in elk geval verlies van alle punten tot dan in de speelbeurt verzameld en bij dubbel 1 gaat tevens de totaalscore terug naar nul. Voor het 1-speler geval, bereken met dynamische programmering:

32

CHAPTER 2. STOCHASTIC DYNAMIC PROGRAMMING • de minimale verwachte waarde van het aantal speelbeurten nodig om 100 punten te bereiken. • Voor N = 7, 10, 15 en 20 de maximale kans om 100 punten te bereiken als de speler over niet meer dan N speelbeurten de beschikking heeft.

10 Een roversbende houdt een bijeenkomst in haar geheime schuilplaats. Buiten staat een agent op wacht die bij toeval de schuilplaats ontdekt heeft en wiens enige doel is de leider van de bende te arresteren. De agent weet dat de schurken uit veiligheidsoverwegingen één voor één naar buiten zullen komen en dat de andere schurken gewaarschuwd zijn en ontsnappen zodra hij één van hen arresteert. De agent wil met maximale kans de leider van de bende die uit 10 leden bestaat arresteren. Gelukkig weet de agent dat de leider de langste van het stel is. Welke strategie moet de agent volgen, wil hij met maximale kans de leider van de bende arresteren? Formuleer dit als een dynamisch-programmeringsprobleem en bereken de optimale strategie. Bepaal de optimale strategie ook voor de situatie dat de bende uit 25, 50 respectievelijk 100 leden bestaat. 11 In de TV show ”De Zwakste Schakel” krijgt een team van spelers een totaal van N vragen voorgelegd, waarbij elke vraag komt uit één van de vraagklasses i = 1, 2, . . . , 9. Na elke correct beantwoorde vraag kan het team besluiten eerst de opgebouwde winst in de huidige keten van correct beantwoorde vragen veilig te stellen alvorens de volgende vraag aan te horen. Zou het team dit doen na een keten van i correct beantwoorde vragen, dan komt in principe een bedrag ci bij het bedrag dat al in de kluis veilig gesteld is en begint een nieuwe keten van vragen. Elke nieuwe keten van vragen begint met een vraag uit vraagklasse 1, is deze vraag goed beantwoord en gaat men door, dan volgt een vraag uit vraagklasse 2, etc. De huisregel is dat het bedrag in de kluis nooit een gegeven topwaarde W mag overschrijden, d.w.z. als reeds een bedrag b in de kluis was en het team halt houdt na een keten van i correct beantwoorde vragen, dan wordt het bedrag in de kluis gelijk aan min(b + ci , W ). Voor de ci ’s geldt 0 < c1 < c2 < · · · < c9 = W , d.w.z een keten van negen goed beantwoorde vragen levert het maximaal te winnen bedrag W op. Zou de huidige keten niet gestopt worden en de volgende vraag wordt foutief beantwoord, dan gaat de in de keten opgebouwde winst verloren en begint een nieuwe keten zolang in totaal nog niet N vragen gesteld zijn. Verondersteld wordt dat voor elke vraag k er een gegeven kans pk is dat het team tot het correcte antwoord op de vraag komt. Bepaal voor het getallenvoorbeeld, c1 = 20, c2 = 50, c3 = 100, c4 = 200, c5 = 300, c6 = 450, c7 = 600, c8 = 800 en c9 = W = 1000 de maximale verwachte opbrengst wanneer N gevarieerd wordt als N = 15, 20 25 en p als p = 0.5, 0.7, 0.75 en 0.8 met pk = p voor alle k. Analyseer ook de structuur van de optimale strategie. Doe het zelfde voor de data c1 = 100, c2 = 250, c3 = 500, c4 = 1000, c5 = 2500, c6 = 5000, c7 = 7500, c8 = 10000 en c9 = W = 12500.

Chapter 3 Appendix 1: Poisson process Een stochastisch proces dat onlosmakelijk verbonden is met de Poisson verdeling is het Poisson proces. Dit is een telproces dat het aantal optredens van een bepaalde gebeurtenis telt over de loop van de tijd. De gebeurtenissen kunnen van velerlei aard zijn: de aankomst van klanten bij een bank, het optreden van ernstige aardbevingen, de aankomst van telefoongesprekken bij een telefooncentrale, de binnenkomst van storingen in een electriciteitscentrale, de binnenkomst van noodoproepen bij een alarmcentrale, etc. In de rest van deze paragraaf zullen we voor het optreden in de tijd van gebeurtenissen de terminologie hanteren van aankomsten van klanten. Wanneer is een telproces een Poisson proces? Daarvoor is het noodzakelijk een onbegrensd grote populatie van potentiële klanten te veronderstellen, waarbij de klanten zich onafhankelijk van elkaar gedragen. Het aankomstproces van klanten bij een bedieningsfaciliteit heet een Poisson proces als het proces de volgende eigenschappen heeft: a. de klanten komen één voor één aan b. de aantallen aankomsten in disjuncte tijdsintervallen zijn onafhankelijk van elkaar c. de kansverdeling van het aantal klanten dat aankomt in een gegeven tijdsinterval heeft een Poisson verdeling waarvan de verwachtingswaarde evenredig is met de lengte van het tijdsinterval.∗ Noteren we met λ = de verwachtingswaarde van het aantal aankomsten in een tijdsinterval van eenheidslengte, ∗

In feite kan de derde aanname verzwakt worden tot de aanname dat de kansverdeling van het aantal aankomsten in een gegeven tijdsinterval alleen afhangt van de lengte van het tijdsinterval en niet van de positie van het interval op de tijdas. Verder moet dan worden aangenomen dat de kans op twee of meer aankomsten in een zeer klein tijdsinterval van lengte h gegeven dat een aankomst in het tijdsinterval heeft plaatsgevonden naar 0 gaat als h naar 0 gaat.

33

34

CHAPTER 3.

APPENDIX 1: POISSON PROCESS

dan stelt eigenschap c. dat voor elke vaste t > 0 geldt P (k aankomsten in een tijdstinterval van lengte t) k −λt (λt) =e voor k = 0, 1, . . . . k! Het getal λ heet de aankomstintensiteit van het Poisson proces. Het Poisson proces geeft een goede modelbeschrijving in veel praktische situaties. De verklaring hiervan is gelegen in het volgende, ruwweg geformuleerde, resultaat. Stel dat op micro niveau sprake is van een zeer groot aantal stochastische processen die onafhankelijk van elkaar zijn en die elk spaarzaam aankomsten genereren. Dan kan aangetoond worden dat op macro niveau de resultante van al deze processen bij goede benadering een Poisson proces is. Dit verklaart bijvoorbeeld waarom het aankomstproces van klanten bij een bank of een postkantoor in veel gevallen goed beschreven wordt door een Poisson proces. Je hebt een zeer grote populatie van potentiële klanten: ieder afzonderlijk persoon gedraagt zich weliswaar volgens een bepaald patroon maar de samenvoeging van al deze typisch onafhankelijke patronen leidt tot een onvoorspelbaar geheel dat zich goed laat beschrijven door een Poisson proces. De eigenschap b. die het ontbreken van na-effecten in het aankomstproces veronderstelt, is de meest kenmerkende eigenschap van het Poisson proces en kan alleen vervuld zijn als de populatie van klanten onbegrensd groot is. Intu¨ıtief zal duidelijk zijn dat deze eigenschap impliceert: voor elk gegeven tijdstip t geldt dat de kansverdeling van de tijd tot de eerstvolgende aankomst niet afhangt van hoe lang geleden het is dat voor het laatst een klant binnenkwam. Deze eigenschap staat bekend als de geheugenloosheidseigenschap van het Poisson proces. Het Poisson proces met deze geheugenloosheidseigenschap staat lijnrecht tegenover een deterministisch aankomst proces waarin klanten met vaste tussentijden aankomen. Voorbeeld A.1 Bij het Centraal Station staan groeptaxi’s. Een groeptaxi vertrekt zodra vier passagiers zijn ingestopt of tien minuten verstreken is sinds de eerste passagier is ingestapt. Veronderstel dat passagiers aankomen volgens een Poisson proces met een gemiddelde van één passagier per drie minuten. (a) Je stapt als eerste passagier in. Wat is de kans dat je tien minuten moet wachten tot het vertrek van de taxi? (b) Je bent als eerste passagier ingestapt en je zit al vijf minuten in de taxi. Inmiddels zijn ook twee andere passagiers ingestapt. Wat is de kans dat je nog vijf minuten moet wachten tot het vertrek van de taxi? Het antwoord op vraag (a) berust op de constatering dat je 10 minuten moet wachten alleen dan als het aantal passagiers dat de komende 10 minuten aankomt

35 minder dan drie is. Kiezen we de minuut als tijdseenheid, dan is dit aantal Poisson verdeeld en heeft de verwachtingswaarde 10λ met λ = 13 . Dus P (je moet 10 minuten wachten) = P (komende 10 minuten stappen 0, 1 of 2 passagiers in) 10λ (10λ)2 = e−10λ + e−10λ + e−10λ = 0.3528. 1! 2! De beantwoording van vraag (b) berust op de geheugenloosheidseigenschap van het Poisson proces. De wachttijd tot de aankomst van de eerstkomende passagier is exponentieel verdeeld met verwachtingswaarde van drie minuten, ongeacht hoe lang geleden de laatste passagier instapte. De kans dat je nog 5 minuten moet wachten is gelijk aan de kans e−5λ = 0.1889 dat de komende vijf minuten er geen nieuwe passagier bij komt. Voorbeeld A.2 Je komt op een random moment tussen 5 uur en kwart over 5 bij de bushalte aan de eerstkomende bus naar huis te nemen. Zowel buslijn 1 als buslijn 3 kun je nemen. Buslijn 1 vertrekt precies op het uur en elke 15 minuten daarna van de halte, terwijl buslijn 3 gemiddeld één keer in de 15 minuten vertrekt waarbij de vertrekmomenten worden bepaald door een Poisson proces. Wat is de kans dat je met buslijn 1 naar huis gaat? Onder de conditie dat vanaf het moment van je aankomst bij de halte het nog x minuten duurt alvorens bus 1 volgens schema aankomt, is de kans dat je met bus 1 naar huis gaat gelijk aan de kans dat de komende x minuten geen bus 3 aankomt. Deze kans is gelijk aan e−λx , ongeacht hoe lang geleden bus 3 voor het 1 laatst aankwam. Hierbij is λ = 15 . Door te middelen over de uniform verdeelde tijd dat het duurt tot bus 1 aankomt, vind je 1 P (je gaat met bus 1 naar huis) = 15

Z

15 0

1 1 e−λx dx = 1 − (> ). e 2

Alternatieve definities van het Poisson proces Bij een Poisson aankomstproces is het aantal aankomsten in een gegeven tijdsinterval discreet verdeeld, maar heeft de tijd tussen twee opeenvolgende aankomsten een continue kansverdeling en wel de exponentiële verdeling.∗ Dit is als ∗

Het Poisson proces heeft zowel een discrete component (Poisson verdeling voor aantal aankomsten) als een continue component(exponentiële verdeling voor tussenaankomsttijden). Deze zaken worden nog wel eens door elkaar gehaald. Een illustratieve situatie is de volgende. Uit fysische experimenten is bekend dat de emissie van alpha-deeltjes uit radioactief materiaal goed beschreven kan worden door een Poisson proces: het aantal deeltjes dat in een gegeven tijdsinterval uitgestoten wordt is Poisson verdeeld en de tijden tussen de uitstoot van opeenvolgende deeltjes zijn onafhankelijk en exponentieel verdeeld.

36

CHAPTER 3.


volgt in te zien: P (tijd tussen twee opeenvolgende aankomsten is groter dan y) = P (in een interval van de lengte y vindt geen aankomst plaats) = e−λy voor elke y > 0. Bij een Poisson aankomstproces met aankomstintensiteit λ heeft de stochastische tijd T tussen twee opeenvolgende aankomsten dus de kansverdelingsfunctie P (T ≤ y) = 1 − e−λy

voor y ≥ 0.

Dit is een continue verdelingsfunctie met als afgeleide de exponentiële kansdichtheidsfunctie λe−λy . De verwachtingswaarde van de exponentiële verdeling is λ1 . Een verwachtingswaarde van λ1 voor de tussenaankomsttijd is in overeenstemming met de betekenis van de aankomstintensiteit λ. Vanwege eigenschap b. van het Poisson proces zal het niemand verbazen dat de tijden tussen de aankomsten van opeenvolgende klanten onafhankelijk van elkaar zijn. De volgende alternatieve definitie van het Poisson proces kan dan worden gegeven: Alternatieve definitie 1. Een aankomstproces van klanten is een Poisson proces als de tussenaankomsttijden van klanten onderling onafhankelijke stochastische variabelen zijn met eenzelfde exponentiële kansverdeling. Deze definitie is de basis voor de simulatie van aankomsttijdstippen in een Poisson aankomstproces met intensiteit λ . Je doet elke keer met de inversie-methode trekkingen uit de exponentiële kansverdelingsfunctie F (t) = 1 − e−λt . Dit doe je door bij een random getal u de tussentijd t op te lossen uit de vergelijking F (t) = 1 − e−λt = u, oftewel t = − λ1 ln(1 − u). Uit de alternatieve definitie van het Poisson proces leiden we vervolgens een derde alternatieve definitie af die voor toepassingen veelal de meest bruikbare is. Daartoe hebben we het begrip ”faalintensiteit” van een positieve stochastische variabele nodig. Dit begrip is het simpelst uit te leggen als we veronderstellen dat de positieve stochast X de levensduur van een machine voorstelt, waarbij we tevens veronderstellen dat X exponentieel verdeeld is met verwachtingswaarde 1 . Dan geldt, ongeacht de waarde van de leeftijd x van de machine µ P (machine van leeftijd x faalt in het komend tijdje ∆x) ≈ µ∆x + o(∆x) for ∆x → 0, waarbij o(∆x) een algemene verzamelnaam is voor een restterm die verwaarloosbaar klein is ten opzichte van ∆x als ∆x zelf heel klein is. Een machine met een exponentiële levensduur heeft dus een constante faalintensiteit µ (“used is good as new”). De afleiding van dit resultaat is als volgt. Op grond van

37 P (X ≤ x) = 1 − e−µx voor x > 0 geldt P (x < X ≤ x + ∆x) P (X > x) −µx e − e−µ(x+∆x) = 1 − e−µ∆x . = 1 − e−µx

P (X ≤ x + ∆x | X > x) =

Uit de reeksontwikkeling 1 − e−u = 1 − u + 1−e

−µ∆x

u2 2!

−

u3 3!

+ · · · volgt

(µ∆x)2 (µ∆x)3 = µ∆x − + − · · · = µ∆x + o(∆x), 2! 3!

waarbij o(∆x) de wiskundige notatie is voor elke functie van ∆x die verwaarloosbaar klein is ten opzichte van ∆x als ∆x zelf heel klein is, dwz o(∆x)/∆x → 0 als ∆x → 0. Op grond van bovenstaande kunnen we nu de derde alternatieve definitie van het Poisson proces geven: Alternative definitie 2. Een aankomstproces {N (t), t ≥ 0} van klanten, waarbij N (t) het totale aantal aankomsten tot tijdstip t aangeeft, heet een Poisson proces met aankomstintensiteit λ als: a. De aantallen aankomsten in disjuncte tijdsintervallen onafhankelijk van elkaar zijn, b. Voor elke t ≥ 0, P (N (t + ∆t) − N (t) = 1) = λ∆t + o(∆t) als ∆t → 0, c. Voor elke t ≥ 0, P (N (t + ∆t) − N (t) ≥ 2) = o(∆t) als ∆t → 0. Deze definitie heeft als bijkomend voordeel dat de definitie eenvoudig uitbreidbaar is naar de situatie waarin de aankomstintensiteit van klanten tijdsafhankelijk is. Samenvoegen en splitsen van Poisson processen In toepassingen komt vaak de situatie voor dat twee Poisson processen samengevoegd worden of dat een Poisson proces uitgedund wordt. Stel dat een call-center telefonische verzoeken voor informatie afhandelt voor twee totaal verschillende bedrijven. Gesprekken voor het ene bedrijf A komen binnen volgens een Poisson proces met aankomstintensiteit λA en onafhankelijk daarvan komen gesprekken voor het andere bedrijf B binnen volgens een Poisson proces met aankomstintensiteit λB . De samenvoeging van deze twee aankomstprocessen is dan een Poisson proces met aankomstintensiteit λA + λB . Verder kan bewezen worden dat een A een gesprek voor bedrijf A is en met random gekozen gesprek met kans λAλ+λ B λB kans λA +λB een gesprek voor bedrijf B. Als voorbeeld voor het splitsen van een Poisson proces beschouw de situatie dat het optreden in de tijd van aardbevingen in een bepaald gebied beschreven

38

CHAPTER 3.


wordt door een Poisson proces met intensiteit λ. Veronderstel dat de sterktes van de aardbevingen onafhankelijk van elkaar zijn, waarbij de sterkte van een aardbeving met kans p geclassificeerd wordt als een zware en met kans 1−p als een niet-zware. Het kan dan bewezen worden dat het optreden in de tijd van zware aardbevingen beschreven wordt door een Poisson proces met intensiteit λp en van niet-zware aardbevingen door een Poisson proces met intensiteit λ(1 − p). Dit resultaat is niet verrassend. Wel verrassend is het resultaat dat de twee Poisson processen die de zware en niet-zware aardbevingen beschrijven onafhankelijk van elkaar zijn. Voorbeeld A.3 Een stuk radioactief materiaal stoot α-deeltjes uit volgens een Poisson proces met een gemiddelde van 0.84 deeltje per seconde. Een Geiger teller registreert, onafhankelijk van elkaar, elk uitgestoten deeltje met kans 0.95. In een 10-seconden periode zijn 12 deeltjes geregistreerd. Wat is de kans dat meer dan 15 deeltjes uitgestoten werden in die periode? Het proces dat de uitstoot van niet-geregistreerde deeltjes beschrijft is een Poisson proces met een intensiteit van 0.84 × 0.05 = 0.0402 deeltje per seconde en dit proces is onafhankelijk van het Poisson proces dat de uitstoot van geregistreerde deeltjes beschrijft. Dus het aantal deeltjes dat door de Geiger teller gemist werd in de 10-seconde periode is Poisson verdeeld met verwachtingswaarde 10 × 0.0402 = 0.402, ongeacht het aantal uitgestote deeltjes dat geregistreerd werd. De gevraagde kans is de kans dat in een gegeven periode van 10 seconde meer dan drie niet-geregistreerde deeltjes worden uitgestoten en is dus gelijk aan P j 1 − 3j=0 e−0.402 (0.402) = 0.00079. j! Voorbeeld A.4 De splitsingseigenschap van het Poisson proces stelt ons in staat om een beroemd resultaat uit de wachttijdtheorie plausibel te maken. Een eachttijdmodel dat vele praktische toepassingen kent, is het zogenoemde M/G/∞ wachttijdmodel met een onbegrensd aantal bediendes en een Poisson aankomstproces van klanten. Dit wachttijdmodel veronderstelt: • klanten komen aan volgens een Poisson proces met een gemiddelde van λ klanten per tijdseenheid, • een onbegrensd aantal bediendes is beschikbaar zodat elke klant bij aankomst direct een vrije bediende krijgt toegewezen, • de bedieningstijd B van een klant heeft een algemene kansverdeling. In het M/G/∞ model is het aantal bezette bediendes dus altijd gelijk aan het aantal aanwezige klanten. De limietverdeling van het aantal bezette bediendes wordt gegeven door j −λE(B) [λE(B)]

lim P (er zijn j bezette bediendes op tijdstip t) = e

t→∞

j!

, j = 0, 1, . . . .

39 In het M/G/∞ model is dus de vorm van de kansverdeling van de bedieningstijd niet van belang voor de limietverdeling van het aantal bezette bediendes, alleen de verwachtingswaarde van de bedieningstijd is relevant. Dit ongevoeligheidsresultaat is uiterst nuttig voor praktische toepassingen. Een illustratie is de volgende. Stel dat olietankers vanuit het Midden-Oosten naar Rotterdam vertrekken volgens een Poisson proces met een gemiddelde van twee per dag. De vaartijden van de tankers worden onafhankelijk van elkaar verondersteld en zijn gamma verdeeld met een verwachting van 12 dagen en een spreiding van anderhalf dag. Het aantal olietankers dat zich op een willekeurig tijdstip op zee bevindt op weg naar Rotterdam is dan Poisson verdeeld met verwachtingswaarde 2 × 12 = 24. De gemiddelde waarde van de vaartijd is alleen relevant en niet de verdeling! Een toepassing van het M/G/∞ model in de voorraadtheorie is als volgt. Voor een exclusief product bij een winkelier wordt het vraagproces beschreven door een Poisson proces met intensiteit λ, d.w.z. klanten komen aan volgens een Poisson proces met intensiteit λ en elke klant vraagt precies één eenheid van het product. Elke keer als een eenheid van het product gevraagd wordt, dan bestelt de winkelier bij de fabriek één eenheid van het product bij. De levertijden van de aanvulorders bij de fabriek zijn onafhankelijk van elkaar en identiek verdeeld met verwachtingswaarde µ. Een vraag van een klant op een moment dat er geen voorraad op de planken is, wordt nageleverd zodra de nabestelling bij de fabriek binnenkomt bij de winkelier. De basisvoorraad van de winkelier is B. Voor dit voorraadmodel geldt dat de som van de voorraad op de planken en het aantal uitstaande aanvulorders bij de fabriek altijd gelijk aan B is. Na enig nadenken volgt dat het M/G∞ model van toepassing is op het aantal uitstaande aanvulorders bij de fabriek zodat de limietverdeling van het aantal uitstaande aanvulorders een Poisson verdeling is met verwachtingswaarde λµ. Hiermee is ook de kans verdeling van de voorraad op de planken bij de winkelier bepaald. Toelichting op ongevoeligheidseigenschap De algemeen geldigheid van bovenstaand resultaat voor de limietverdeling van het aantal bezette bediendes kunnen we plausibel maken door een probabilistisch bewijs te geven voor het geval de bedieningstijd slechts een eindig aantal waarden kan aannemen. Beschouw eerst het speciale geval van een constante bedieningstijd D = E(B). Het resultaat volgt dan direct door op te merken dat voor elk tijdstip t met t > D geldt dat het aantal bezette bediendes op tijdstip t gelijk moet zijn aan het aantal klanten dat in het tijdsinterval (t − D, t] is aangekomen (immers elke klant al in bediening op tijdstip t − D heeft het systeem op tijdstip t verlaten en elke klant die tussen t − D en t binnenkomt is nog in bediening op tijdstip t). Het aantal klanten dat in een tijdsinterval (t − D, t] van de lengte D aankomt, heeft een Poisson verdeling met verwachtingswaarde λD, . Dus voor het het geval van een deterministische bedieningstijd D geldt voor elk tijdstip t > D dat de kansverdeling van het aantal bezette

40

CHAPTER 3.


bediendes Poisson verdeeld is met verwachtingswaarde λD. Deze afleiding laat zich direct uitbreiden tot de situatie dat de bedieningstijd S een eindig aantal waarden D1 , . . . , Dm kan aannemen met respectievelijke kansen p1 , . . . , pm . Voor het gemak nemen we m = 2. Klanten met bedieningstijd D1 noemen we klanten van type 1 en klanten met bedieningstijd D2 noemen we klanten van type 2. Op grond van een splitsingseigenschap van het Poisson proces komen klanten van types 1 en types 2 aan volgens onafhankelijke Poisson processen met respectievelijke parameters λ1 = λp1 en λ2 = λp2 . Kies nu t vast zodat t > D1 en t > D2 . Bovenstaande redenering geeft dat op tijdstip t het aantal klanten van type 1 in bediening Poisson verdeeld is met verwachtingswaarde λ1 D1 en het aantal klanten van type 2 in bediening Poisson verdeeld is met verwachtingswaarde λp2 . De som van twee onafhankelijke Poisson verdeelde stochasten is weer Poisson verdeeld. Dus het aantal bezette bedienden op tijdstip t is Poisson verdeeld met verwachtingswaarde λ1 D1 + λ2 D2 = λ(p1 D1 + p2 D2 ) = λE(S). Niet-homogeen Poisson proces In veel praktische situaties is de aankomstintensiteit van klanten tijdsafhankelijk. De alternatieve definitie van het Poisson proces, waarin aankomsten worden geteld door zeer kleine tijdsintervallen te beschouwen, laat zich direct uitbreiden naar het geval van tijdsafhankelijke aankomsten. Een telproces {N (t), t ≥ 0}, waarbij N (t) het totale aantal aankomsten tot tijdstip t aangeeft, heet een niethomogeen Poisson proces met aankomstintensiteitsfunctie λ(t) als: a. De aantallen aankomsten in disjuncte tijdsintervallen onafhankelijk van elkaar zijn, b. Voor elke t ≥ 0, P (N (t + ∆t) − N (t) = 1) = λ(t)∆t + o(∆t) als ∆t → 0, c. Voor elke t ≥ 0, P (N (t + ∆t) − N (t) ≥ 2) = o(∆t) als ∆t → 0. Hierbij is o(h) een generieke notatie voor elke functie f (h) met de eigenschap dat f (h) sneller naar nul gaat dan h zelf als h naar nul gaat (bijvoorbeeld f (h) = h2 is een o(h) functie). Preciezer gesteld, een functie f (h) is een o(h) functie voor h → 0 als limh→0 f (h)/h = 0. Voor een niet-homogeen Poisson aankomstproces {N ((t), t ≥ 0} kan worden aangetoond dat voor alle t1 , t2 ≥ 0 met t2 > t1 geldt dat het totale aantal aankomsten in het tijdsinterval (t1 , t2 ) een Poisson verdeling heeft met verwachtingswaarde M (t2 ) − M (t1 ), waarbij Z t M (t) = λ(x)dx voor t ≥ 0. 0

De functie M (t) geeft dus het verwachte aantal aankomsten in (0, t). Voorbeeld A.5 Bij een pinautomaat komen klanten aan volgens een niethomogeen Poisson proces. Nemen we zes uur ’s ochtends als t = 0 en het uur als

41 tijdseenheid dan komen klanten tussen zes uur ’s ochtends en 12 uur ’s middags aan met intensiteit λ(t) = 7.5t voor 0 ≤ t ≤ 6 en tussen 12 uur ’s middags en zes uur de volgende ochtend met intensiteit λ(t) = 45 − 2.5(t − 6) voor 6 ≤ t ≤ 24. Wat is de kansverdeling van het totale aantal klanten dat per etmaal van de pinautomaat gebruik maakt? Het antwoord is dat dit aantal een Poisson verdeling heeft met verwachtingswaarde Z 24 Z 6 Z 24 λ(t)dt = 7.5tdt + [45 − 2.5(t − 6)]dt = 540. 0

0

6

Stel dat elke klant 50 euro pint met kans 34 en 100 euro met kans 41 . Wat is de verwachting van het totale bedrag dat per etmaal wordt gepind? Noem een klant voor 50 euro een type 1 klant en een klant voor 100 euro een type 2 klant. De type 1 klanten komen dan aan volgens een niet-homogeen Poisson proces met intensiteit 34 λ(t) en de type 2 klanten volgens een niet-homogeen Poisson proces met intensiteit 14 λ(t) (waarom?). Het totale aantal keer dat 50 euro per etmaal gepind wordt is dan Poisson verdeeld met verwachtingswaarde 34 × 540 = 405 en het totale aantal keer dat 100 euro ge¨ınd wordt is Poisson verdeeld met verwachtingswaarde 41 × 540 = 135. Het totale bedrag dat per etmaal gepind wordt, heeft dus verwachtingswaarde µ = 50 × 405 + 100 × 135 = 33 750 euro. Voor het geval de aankomstintensiteitsfunctie λ(t) begrensd in t ≥ 0 is, kan een interessante alternatieve beschrijving van het niet-homogeen Poisson aankomstproces gegeven worden. Kies een vast getal λ met λ ≥ λ(t) voor alle t ≥ 0. Construeer nu als volgt een telproces: 1. Aankomsten vinden plaats volgens een Poisson proces met intensiteit λ. 2. Een aankomst die plaatsvindt op tijdstip s wordt geaccepteerd met kans λ(s)/λ en wordt verworpen anders. Het telproces {N (t), t ≥ 0} met N (t) het aantal geaccepteerde aankomsten tot tijdstip t is dan een niet-homogeen Poisson proces met aankomstintensiteitsfunctie λ(t). Dit is als volgt in te zien. De kans dat het Poisson proces een aankomst genereert in (t, t + ∆t) is λ∆t + o(∆t) voor ∆t klein en de kans dat deze aankomst geaccepteerd wordt is dan λ(t)/λ, oftewel de kans dat in (t, t+∆t) een geaccepteerde aankomst plaatsvindt is gelijk aan (λ∆t) × (λ(t)/λ) + o(∆t) = λ(t)∆t + o(∆t). Deze constructie van een niet-homogeen Poisson proces is nuttig voor het simuleren van aankomsten in een niet-homogeen Poisson proces.

42

CHAPTER 3.


Chapter 4 Appendix 2: Renewal-reward processes 4.1

Basic theory

A powerful tool in the analysis of numerous applied probability models is the renewal-reward model. This model is also very useful for theoretical purposes. Limit theorems for Markov chains can be proved by using the renewal-reward theorem. The renewal-reward model is a simple and intuitively appealing model that deals with a so-called regenerative process on which a cost or reward structure is imposed. Many stochastic processes have the property of regenerating themselves at certain points in time so that the behaviour of the process after the regeneration epoch is a probabilistic replica of the behaviour starting at time zero and is independent of the behaviour before the regeneration epoch. A formal definition of a regenerative process is as follows. Definition 2.2.1 A stochastic process {X(t), t ∈ T } with time-index set T is said to be regenerative if there exists a (random) epoch S1 such that: (a) {X(t + S1 ), t ∈ T } is independent of {X(t), 0 ≤ t < S1 }, (b) {X(t + S1 ), t ∈ T } has the same distribution as {X(t), t ∈ T }. It is assumed that the index set T is either the interval T = [0, ∞) or the countable set T = {0, 1, ...}. In the former case we have a continuous-time regenerative process and in the other case a discrete-time regenerative process. The state space of the process {X(t)} is assumed to be a subset of some Euclidean space. The existence of the regeneration epoch S1 implies the existence of further regeneration epochs S2 , S3 , ... having the same property as S1 . Intuitively speaking, a regenerative process can be split into independent and identically distributed renewal cycles. A cycle is defined as the time interval between two consecutive regeneration epochs. Examples of regenerative processes are: (i) The continuous-time process {X(t), t ≥ 0} with X(t) denoting the number 43

44

CHAPTER 4.

APPENDIX 2: RENEWAL-REWARD PROCESSES

of customers present at time t in a single-server queue in which the customers arrive according to a renewal process and the service times are independent and identically distributed random variables. It is assumed that at epoch 0 a customer arrives at an empty system. The regeneration epochs S1 , S2 , ... are the epochs at which an arriving customer finds the system empty. (ii) The discrete-time process {In , n = 0, 1, ...} is a Markov chain with a finite state space and a state 0 that can be reached from any other state. The regeneration epochs are the epochs at which the Markov chain makes a transition to state 0. Let us define the random variables Cn = Sn − Sn−1 , n = 1, 2, ..., where S0 = 0 by convention. The random variables C1 , C2 , ... are independent and identically distributed. In fact the sequence {C1 , C2 , ...} underlies a so-called renewal process in which the events are the occurrences of the regeneration epochs. The random variable Cn can be interpreted as Cn = the length of the nth renewal cycle,

n = 1, 2, ... .

Note that the cycle length Cn assumes values from the index set T. In the following it is assumed that 0 < E(C1 ) < ∞. In many practical situations a reward structure is imposed on the regenerative process {X(t), t ∈ T }. The reward structure usually consists of reward rates that are earned continuously over time and lump rewards that are only earned at certain state transitions. Denote by Rn = the total reward earned in the nth renewal cycle,

n = 1, 2, ... .

It is assumed that R1 , R2 , ... are independent and identically distributed random variables. In applications Rn typically depends on Cn . In case Rn can take on both positive and negative values, it is assumed that E(|R1 |) < ∞. Let R(t) = the cumulative reward earned up to time t. The process {R(t), t ≥ 0} is called a renewal-reward process. We are now ready to prove a theorem of utmost importance. Theorem 2.2.1 (renewal-reward theorem) E(R1 ) R(t) = t→∞ t E(C1 ) lim

with probability 1.

In other words, for almost any realization of the process, the long-run average reward per time unit is equal to the expected reward earned during one cycle divided by the expected length of one cycle.

4.1. BASIC THEORY

45

To prove this theorem we first establish the following lemma: Lemma 2.2.1 For any t ≥ 0, let N (t) be the number of cycles completed up to time t. Then N (t) 1 lim = with probability 1. t→∞ t E(C1 ) Proof By the definition of N (t), we have C1 + ... + CN (t) ≤ t < C1 + ... + CN (t)+1 . Since P {C1 + ... + Cn < ∞} = 1 for all n ≥ 1, it is not difficult to verify that lim N (t) = ∞

t→∞

with probability 1.

The above inequality gives C1 + ... + CN (t) C1 + ... + CN (t)+1 N (t) + 1 t ≤ < . N (t) N (t) N (t) + 1 N (t) By the strong law of large numbers for a sequence of independent and identically distributed random variables, we have C1 + ... + Cn = E(C1 ) n→∞ n lim

with probability 1.

Hence, by letting t → ∞ in the above inequality, the desired result follows. Lemma 2.2.1 is also valid when E(C1 ) = ∞ provided that P {C1 < ∞} = 1. The reason is that the strong law of large numbers for a sequence {Cn } of nonnegative random variables does not require that E(C1 ) < ∞. Next we prove Theorem 2.2.1. Proof of Theorem 2.2.1 For ease, let us first assume that the rewards are non-negative. Then, for any t > 0, N (t) X

N (t)+1

Ri ≤ R(t) ≤

i=1

X

Ri .

i=1

This gives PN (t)

PN (t)+1 Ri N (t) Ri N (t) + 1 R(t) × ≤ ≤ i=1 × . N (t) t t N (t) + 1 t i=1

By the strong law of large numbers for the sequence {Rn }, we have n

1X Ri = E(R1 ) n→∞ n i=1 lim

with probability 1.

46

CHAPTER 4.


As pointed out in the proof of Lemma 2.2.1, N (t) → ∞ with probability 1 as t → ∞. Letting t → ∞ in the above inequality and using Lemma 2.2.2 the desired result next follows for the case that the rewards are non-negative. In case the rewards can assume both positive and negative values, then the theorem is proved by treating the positive and negative parts of the rewards separately. We omit the details. In a natural way Theorem 2.2.1 relates the behaviour of the renewal-reward process over time to the behaviour of the process over a single renewal cycle. It is noteworthy that the outcome of the long-run average actual reward per time unit can be predicted with probability 1. If we are going to run the process over an infinitely long period of time, then we can say beforehand that in the long run the average actual reward per time unit will be equal to the constant E(R1 )/E(C1 ) with probability 1. This is a much stronger and more useful statement than the statement that the long-run expected average reward per time unit equals E(R1 )/E(C1 ) (it indeed holds that limt→∞ E[R(t)]/t = E(R1 )/E(C1 ); this expected-value version of the renewal-reward theorem is a direct consequence of Theorem 2.2.1 when R(t)/t is bounded in t but otherwise requires a hard proof). Also it is noted that for the case of non-negative rewards Rn the renewal-reward theorem is also valid when E(R1 ) = ∞ (the assumption E(C1 ) < ∞ cannot be dropped for Theorem 2.2.1). Example 2.2.1. Alternating up- and downtimes Suppose a machine is alternately ’up and down’. Denote by U1 , U2 , ... the lengths of the successive up-periods and by D1 , D2 , ... the lengths of the successive down-periods. It is assumed that both {Un } and {Dn } are sequences of independent and identically distributed random variables with finite positive expectations. The sequences {Un } and {Dn } are not required to be independent of each other. Assume that an up-period starts at epoch 0. What is the long-run fraction of time the machine is down? The answer is the long-run fraction of time the machine is down E(D1 ) = with probability 1. E(U1 ) + E(D1 ) It is convenient to define first a regenerative process that is appropriate for the situation considered. For any t > 0, let the random variable ½ 1 if the machine is up at time t X(t) = 0 otherwise. The continuous-time process {X(t)} is a regenerative process.∗ The epochs at which an up-period starts can be taken as regeneration epochs. The long-run ∗

In the state description of the process {X(t)} you could also have included information about the elapsed up-time or the elapsed down-time. The expanded process would then have possessed the Markov property, but this is not relevant for having a regenerative process.

4.1. BASIC THEORY

47

fraction of time the machine is down can be interpreted as a long-run average cost per time unit by assuming that a cost at rate 1 is incurred while the machine is down and a cost at rate 0 otherwise. A regeneration cycle consists of an upperiod and a down-period. Hence E(length of one cycle) = E(U1 + D1 ) and E(cost incurred during one cycle) = E(D1 ). By applying the renewal-reward theorem, it follows that the long-run average cost per time unit equals E(D1 )/[E(U1 ) + E(D1 )] proving the result. Remark. In case the up-time or down-time has a (continuous) probability density, then E(D1 ) lim P {the machine is up at time tB} = . t→∞ E(U1 ) + E(D1 ) The intermediate step of interpreting the long-run fraction of time that the process is in a certain state as a long-run average cost (reward) per time unit is very helpful in many situations. Limit theorems for regenerative processes An important application of the renewal-reward theorem is the characterization of the long-run fraction of time a regenerative process {X(t), t ∈ T } spends in some given set B of states. For the set B of states, define for any t ∈ T the indicator variable ½ 1 if X(t) ∈ B, IB (t) = 0 if X(t) ∈ / B. Also, define the random variable TB = the amount of time the process spends in the set B of states during one cycle. RS Note that TB = 0 1 IB (u)du for a continuous-time process {X(t)}; otherwise, TB equals the number of indices 0 ≤ k ≤ S1 with X(k) ∈ B. The following theorem is an immediate consequence of the renewal-reward theorem. Theorem 2.2.3 For the regenerative process {X(t)} holds: the long-run fraction of time the process spends in the set B of states E(TB ) = with probability 1. E(C1 ) That is, 1 lim t→∞ t

Z

t

IB (u)du = 0

E(TB ) E(C1 )

with probability 1

48

CHAPTER 4.


for a continuous-time process {X(t)} and n

E(TB ) 1X IB (k) = n→∞ n E(C1 ) k=0 lim

with probability 1

for a discrete-time process {X(n)}. Proof The long-run fraction of time the process {X(t)} spends in the set B of states can be interpreted as a long-run average reward per time unit by assuming that a reward at rate 1 is earned while the process is in the set B and a reward at rate 0 is earned otherwise. Then E(reward earned during one cycle) = E(TB ). The desired result next follows by applying the renewal-reward theorem. Since E(IB (t)) = P {X(t) ∈ B}, we have as consequence of Theorem 2.2.3 and the bounded convergence theorem that Z E(TB ) 1 t P {X (u) ∈ B}du = . lim t→∞ t 0 E(C1 ) Rt Note that (1/t) 0 P {X (u) ∈ B}du can be interpreted as the probability that an outside observer arriving at a randomly chosen point in (0, t) finds the process in the set B. Remark. In case an aperiodicity condition is satisfied (e.g., the length of a cycle has a continuous probability density on some interval), we also have that lim P {X (t) ∈ B} =

t→∞

E(TB ) . E(C1 )

This ordinary limit need not always exist. A counterexample is provided by periodic discrete-time Markov chains.

4.2

Poisson arrivals see time averages

In the analysis of queueing (and other) problems one sometimes needs the longrun fraction of time the system is in a given state and sometimes needs the long-run fraction of arrivals who find the system in a given state. These averages can often be related to each other, but in general they are not equal to each other. To illustrate that the two averages are in general not equal to each other, suppose that customers arrive at a service facility according to a deterministic process in which the interarrival times are one minute. If the service of each customer is uniformly distributed between 14 minute and 34 minute, then the long-run fraction of time the system is empty equals 12 , whereas the long-run fraction of arrivals

4.2. POISSON ARRIVALS SEE TIME AVERAGES

49

finding the system empty equals 1. However the two averages would have been the same if the arrival process of customers had been a Poisson process. As a prelude to the generally valid property that Poisson arrivals see time averages, we first analyze two specific problems by the renewal-reward theorem. Example 2.2.2. A manufacturing problem Suppose that jobs arrive at a work station according to a Poisson process with rate λ. The work station has no buffer to store temporarily arriving jobs. An arriving job is accepted only when the work station is idle, and is lost otherwise. The processing times of the jobs are independent random variables having a common probability distribution with finite mean β. What is the long-run fraction of time the work station is busy and what is the long-run fraction of jobs that are lost? These two questions are easily answered by using the renewal-reward theorem. Let us define the following random variables. For any t ≥ 0, let ½ 1 if the work station is busy at time t. I(t) = 0 otherwise. Also, for any n = 1, 2, ..., let ½ 1 if the work station is busy just prior to the nth arrival In = 0 otherwise. The continuous-time process {I(t)} and the discrete-time process {In } are both regenerative. The arrival epochs occurring when the working station is idle are regeneration epochs for the two processes (why?). Let us say that a cycle starts each time an arriving job finds the working station idle. The long-run fraction of time the working station is busy is equal to the expected amount of time the working station is busy during one cycle divided by the expected length of one cycle. The expected length of this busy period equals β. Since the Poisson arrival process is memoryless, the expected length of the idle period during one cycle equals the mean interarrival time 1/λ. Hence, with probability 1, the long-run fraction of time the work station is busy β . = β + 1/λ The long-run fraction of jobs that are lost equals the expected number of jobs lost during one cycle divided by the expected number of jobs arriving during one cycle. Since the arrival process is a Poisson process, the expected number of (lost) arrivals during the busy periods equals λ × E(processing time of a job) = λβ. Hence, with probability 1, the long-run fraction of jobs that are lost λβ . = 1 + λβ

50

CHAPTER 4.


Thus, we obtain from the remarkable result the long-run fraction of arrivals finding the work station busy = the long-run fraction of time the work station is busy . Incidentally, it is interesting to note that in this loss system the long-run fraction of lost jobs is insensitive to the form of the distribution function of the processing time but needs only the first moment of this distribution. This simple loss system is a special case of the famous Erlang loss model. Example 2.2.3. An inventory model Consider a single-product inventory system in which customers asking for the product arrive according to a Poisson process with rate λ. Each customer asks for one unit of the product. Each demand which cannot be satisfied directly from stock on hand is lost. Opportunities to replenish the inventory occur according to a Poisson process with rate µ. This process is assumed to be independent of the demand process. For technical reasons a replenishment can only be made when the inventory is zero. The inventory on hand is raised to the level Q each time a replenishment is done. What is the long-run fraction of time the system is out of stock? What is the long-run fraction of demand that is lost? In the same way as in Example 2.4.1, we define the random variables ½ 1 if the system is out of stock at time t I(t) = 0 otherwise and In =

½

1 0

if the system is out of stock when the nth demand occurs otherwise.

The continuous-time process {I(t)} and the discrete-time process {In } are both regenerative. The regeneration epochs are the demand epochs at which the stock on hand drops to zero (why?). Let us say that a cycle starts each time the stock on hand drops to zero. The system is out of stock during the time elapsed from the beginning of a cycle until the next inventory replenishment. This amount of time is exponentially distributed with mean 1/µ. The expected amount of time it takes to go from stock level Q to 0 equals Q/λ. Hence, with probability 1, the long-run fraction of time the system is out of stock 1/µ . = 1/µ + Q/λ To find the fraction of demand that is lost, note that the expected amount of demand lost in one cycle equals λ × E(amount of time the system is out of stock during one cycle) = λ/µ. Hence, with probability 1, the long-run fraction of demand that is lost λ/µ . = λ/µ + Q

4.3. PROBLEMS

51

Again we find the remarkable result: the long-run fraction of customers finding the system out of stock = the long-run fraction of time the system is out of stock. The remarkable identies we have found in the two examples are particular instances of the property ‘Poisson arrivals see time averages’. Roughly stated, this property expresses that in statistical equilibrium the distribution of the state of the system just prior to an arrival epoch is the same as the distribution of the state of the system at an arbitrary epoch when arrivals occur according to a Poisson process. An intuitive explanation of the property ‘Poisson arrivals see time averages’ is that Poisson arrivals occur completely randomly in time.

4.3

Problems

1. The lifetime of a street lamp is a continuous random variable L with probability distribution function F (x) and probability density f (x). The street lamp is replaced by a new one upon failure or upon reaching the critical age a for a given a, whichever occurs first. A cost of cf > 0 is incurred for each failure replacement and a cost of cp > 0 for each preventive replacement, where cp < cf . The lifetimes of the street lamps are independent of each other. • Define a regenerative process and specify its regeneration epochs. • Show that the long-run average cost per time R a unit under the age-replacement rule equals g(a) = [cp + (cf − cp )F (a)] / 0 {1 − F (x)}dx. 2. At a production facility orders arrive according to a renewal process with a mean interarrival time 1/λ. A production is started only when N orders have accumulated. The production time is negligible. A fixed cost of K > 0 is incurred for each production setup and holding costs are incurred at the rate of hj when j orders are waiting to be processed • Define a regenerative stochastic process and identify its regeneration epochs. • Determine the long-run average cost per unit time. • What value of N minimizes the long-run average cost per unit time? 3. Passengers arrive at a bus stop according to a Poisson process with rate λ. Busses depart from the stop according to a renewal process with interdeparture time A. Using renewal-reward processes, prove that the long-run average waiting time per passenger equals E(A2 )/2E(A). (Hint: imagine a cost structure for which the system incurs a cost at rate x when x is the time until the next 2 and is a bus arrival). The coefficient of variation c2A is defined by c2A = Eσ 2(A) (A)

52

CHAPTER 4.


dimensionless quantity. Verify that the average waiting time E(A2 )/2E(A) per passenger can be written as 21 (1 + c2A )E(A). Thus the average waiting time of a passenger is larger than the average interdeparture time of busses if c2A > 1, that is, if busses arrive irregularly. This is known as the waiting-time paradox. 4. A common car-service between cities in Israel is called ’Sheroot’. A Sheroot is a 7-seat cab that leaves from its stand as soon as it has collected seven passengers. Suppose that potential passengers arrive at the stand according to a Poisson process with rate λ. An arriving person who sees no cab at the stand goes elsewhere and is lost for the particular car-service. Empty cabs pass the stand according to a Poisson process with rate µ. An empty cab stops only at the stand when there is no other cab. • Define a regenerative process and identify its regeneration epochs. • Determine the long-run fraction of time there is no cab at the stand and determine the long-run fraction of customers who are lost. Can you intuitively explain why these two fractions are equal to each other? 5. Big Jim, a man of few words, runs a one-man business. This business is called upon by loansharks to collect overdue loans. Big Jim takes his profession seriously and accepts only one assignment at a time. Assignments that are refused by Jim when het is working on another assignment are handled by a colleague of Big Jim. The assignments are classified by Jim into n different categories j = 1, ..., n. An assignment of type j takes him a random number of τj days and gives a random profit of ξj dollars for j = 1, ..., n. Assignments of the types 1, ..., n arrive according to independent Poisson processes with respective rates λ1 , ..., λn . • Define a regenerative process and identify its regeneration epochs. • Determine the long-run average pay-off per time unit for Big Jim. • Determine the long-run fraction of time Big Jim is at work and the longrun fraction of the assignments that are not accepted. Can you intuitively explain why these two fractions are equal to each other?

Aanvullend dictaat Stochastische Operations Research I. H.C. Tijms

Recommend Documents