Betrouwbaarheid Een simulatie beoogt één of i.h.a. twee of meerdere scenario’s te evalueren en te vergelijken, bij •
Monte Carlo (MC) simulatie voor een groot aantal instelwaarden, voor één of meerdere parameters. Hierbij stelt elke instelling van waarden één scenario voor.
•
Logistieke Simulaties op basis van discrete event (DE) simulatie voor een relatief klein aantal, bijvoorbeeld slechts twee, scenario’s zoals met verschillende inrichtingen of volgordelijkheden van werkzaamheden of capaciteiten.
In beide gevallen (MC of DE) levert simulatie in principe voor één scenario één getal op voor elk gekozen prestatiemaat, bijvoorbeeld een doorlooptijd of de kosten. Dus alsof men voor dat scenario één keer met een dobbelsteen heeft gegooid. Maar wat zegt dit, een 1, 2, 3, 4, 5, of 6 als uitkomst? De voor de hand liggende aanpak is derhalve, deze dobbelsteen (dit scenario) opnieuw te gooien (te simuleren) en na een redelijk groot aantal herhalingen de uitkomst te middelen. Het hieruit bepaalde getal, bijvoorbeeld 3.71 na 30 keer dobbelen, kan men dan terecht wel als redelijk representatieve uitkomst zien. Bij simulatie van een scenario wordt een dergelijke uitkomst door “middeling” in de regel als dé uitkomst voor het gesimuleerde scenario gezien. Hierbij moet men verschillende scenario’s zien als verschillende niet-identieke dobbelstenen, zeg gekleurde dobbelstenen met verschillende kansen. Maar wat is redelijk? Hoeveel waarde mag men aan één zo’n getal hechten? Of preciezer: hoeveel betrouwbaarheid mag men hieraan toekennen? Dit is bij simulatie studies in de praktijk een veelal vergeten en onbeantwoorde vraag, alvorens men daadwerkelijk een simulatie en een vergelijk van scenario’s als representatief en betrouwbaar mag zien.
Betrouwbaarheidsinterval Standaard statistiek lijkt hierop een passend antwoord te geven. De wijze waarop wordt hieronder eerst beknopt uiteengezet. Laat de kansvariabelen (stochastische grootheden) X1, ..., Xn een n-tal ónafhankelijke waarnemingen (herhalingen) voorstellen van één en dezelfde kansvariabele (bijvoorbeeld het aantal ogen van een dobbelsteen) met verwachtingswaarde µ en standaarddeviatie σ. Voor n voldoende groot, praktisch in de orde van n ≥ 50, geldt dan dat het zogenaamde steekproefgemiddelde
X=
X1 + K + X n n
bij redelijke benadering normaal verdeeld is met eveneens verwachting µ maar met standaarddeviatie σ / n . (Dit berust op de zogenaamde centrale limietstelling.) Op grond hiervan geldt dat bij daadwerkelijke individuele waarnemingen (realisaties):
X 1 = x1 , X 2 = x2 , …, X n = xn de gemiddelde waarde x = ( x1 + K + xn ) / n
met kans (1-α) 100% bevat is in het interval:
σ σ , µ + z(1−α / 2) µ − z(1−α / 2) . n n Hierbij is z(1-α/2) het welbekende z-punt zodanig dat voor de standaard normale verdeling Z geldt (zie plaatje voor α=.67 en α=.95)
P ( − z(1−α / 2) < Z < z(1−α / 2) ) = (1 − α ) .
Voor (1-α) 100% = 95% bijvoorbeeld geldt z(1-α/2) = 1.96. Wanneer σ niet bekend is kan men deze schatten en vervangen door de waarde s=
1 2 2 ( x1 − x ) + ... + ( xn − x ) n −1
en dient het z-punt vervangen te worden door het corresponderende tn-1,(1-α/2)-punt voor een zogenaamde Student verdeling met (n-1) vrijheidsgraden (feitelijk gebaseerd op een onderliggende aanname van normaliteit). Voor n redelijk groot, zeg n ≥ 50, en (1-α)100% = 95% is het gerechtvaardigd het z- of t-punt simpelweg af te schatten met de waarde 2. Oftewel, bij een geschatte waarde S = s geldt dat met 95% betrouwbaarheid: s s x ∈ µ − 2 , µ +2 . n n Dus in woorden, dat met 95% betrouwbaarheid:
x niet meer dan 2
s af ligt van µ, n
oftwel, puur taalkundig vertaald, dat met 95% betrouwbaarheid:
µ niet meer dan 2
s af ligt van x . n
Dus, wanneer µ onbekend is kan op basis van een n-voudige steekproef (of simulatie) met waarnemingen x1, ..., xn een 95%-betrouwbaarheidsinterval opgesteld worden als:
µ ∈ x − 2
s s , x +2 . n n
(1)
Let op (Interpretatie) Dit is echter een gevaarlijke, zo niet feitelijk onjuiste uitspraak en moet zorgvuldig geïnterpreteerd worden. Deze formulering suggereert namelijk alsof de verwachtingswaarde µ (zoals 3.5 voor een zuivere dobbelsteen) aan onzekerheid onderhevig is. Dit is geenszins het geval. Integendeel, µ representeert een vast, zij het onbekend, getal waarvoor men o.b.v. de steekproef (of simulatie) een schatting wil geven.
De juiste interpretatie van (1) behoort te zijn: het aldus opgestelde (betrouwbaarheids) interval in (1) o.b.v. het steekproefgemiddelde x (en de steekproef standaarddeviatie s) zal in 95% van de gevallen de échte verwachtingswaarde µ bevatten. Anders gesteld: bij herhaalde uitvoering van een dergelijke steekproef, elk met n waarnemingen, bevat in 95 % van de gevallen het bij die steekproef behorende betrouwbaarheidsinterval zoals boven opgesteld, de echte waarde µ (zie onderstaand figuur). [
[
*µ
]
]
Experiment 1 144444 244444 3 Experiment32 144444 244444 95%-BI B1
95%-BI B2
Ondanks deze oneigenlijke (mis)interpretatie wordt i.h.a. toch simpelweg gesproken van een betrouwbaarheidsinterval voor µ. Ook wel kortweg aangeduid met: 95%-BI. Door n voldoende groot te kiezen (waarbij vooraf een pilot experiment kan worden uitgevoerd om een eerste indicatie van zowel x als s te verkrijgen) kan aldus een gewenste scherpte voor het betrouwbaarheidsinterval worden verkregen. Op basis van deze (95 %) betrouwbaarheidsintervallen kan nu ook aan simulatieresultaten een mate van betrouwbaarheid worden toegekend. Daarbij dient een belangrijk onderscheid gemaakt te worden tussen •
MC-simulatie
•
Logistieke (DE) Simulatie
I Monte Carlo simulatie
In dit geval is toepassing van bovenstaande betrouwbaarheidsinterval direct en 100% gerechtvaardigd. Voor elk scenario, zeg van in totaal N scenario’s, van instelwaarden, voert men n keer één simulatie uit (telkens door het opnieuw trekken uit de diverse benodigde verdelingen). Op basis hiervan stelt men voor dit scenario een 95%-BI op voor de gewenste prestatie (bijvoorbeeld opbrengst). Omdat de trekkingen zowel binnen elk van de scenario’s als voor de verschillende scenario’s onderling onafhankelijk dienen te worden genomen, is aan de onderliggende aanname van onafhankelijkheid voldaan en zijn de BI gerechtvaardigd. In onderstaand figuur is dit geschetst in het geval van één instelwaarde, waarbij B1, B2, …, Bn de 95%-BI voorstellen.
-
B2 B1
-
-
-
- -
-
-
BN
-
1 2
N
Instelwaarde
II Logistieke (Discrete Event) Simulatie
Voor dit aanmerkelijk vaker voorkomende type van computer simulatie ligt het toepassen van de standaard statistisch gerechtvaardigde betrouwbaarheidsintervallen gecompliceerder. Er is namelijk bij een in de tijd doorlopende simulatie geen sprake van duidelijk ónafhankelijke waarnemingen. Als bijvoorbeeld in een wachtrij simulatie de wachttijd van klantnummer 943 groot is, mag men dit ook verwachten voor klantnummer 942 of 944. Voor het opstellen van een BI echter is de onafhankelijkheid van waarnemingen wezenlijk. Een specifiekere opzet voor rechtvaardiging van ónafhankelijke ‘waarnemingen’ is dus vereist. Verschillende methoden zijn hiertoe voor handen. Deze worden hieronder kort uiteengezet aan de hand van een wachtrij simulatie. In de eenvoudigste situatie komen klanten binnen en sluiten aan bij de rij voor het loket (of gaan direct naar het loket als deze vrij is). We zijn geïnteresseerd in de verwachte verblijftijd van een klant die op een willekeurig moment binnenkomt. Bij het verlaten wordt daartoe de verblijftijd van de klant gemeten. Duidelijk is dat de verblijftijd van een klant sterk gecorreleerd is met zijn voorganger(s) en de klanten die na hem komen. De hierboven beschreven methode kan dus niet worden gebruikt. 1.
Replicatie methode
Bij de replicatie methode wordt een simulatie gestart en een vast aantal waardes van de grootheid waarin we zijn geïnteresseerd geregistreerd. Hiervan wordt het gemiddelde bepaald. Bijvoorbeeld voor een wachtrij simulatie wordt de gemiddelde wachttijd van 10000 klanten berekend. Dit levert ons de eerste waarneming x1 op. Dit proces met 10000 klanten wordt n keer herhaald wat leidt n waarnemingen x1, …, xn. Als nu elke simulatie run (van 10000 klanten) met dezelfde startsituatie begint dan zijn deze waarnemingen onafhankelijk van elkaar en identiek verdeeld. Het eenvoudigste is het om de
simulatie te starten met een leeg systeem, oftewel de eerste klant hoeft nooit te wachten. Dit heeft bovendien een positieve invloed op de verblijftijd van de volgende klanten. De waarnemingen x1, …, xn worden dus beïnvloed door de gekozen startsituatie. Dit wordt de ‘initial bias’ genoemd. Om dit effect te vermijden worden de waarden tijdens het eerste gedeelte van de simulatierun uitgesloten in de berekeningen. Dit gedeelte dient slechts ter “opwarming” van de simulatie. Bijvoorbeeld in de wachtrij simulatie wordt de gemiddelde verblijftijd bepaald aaan de hand van de verblijftijden van de klanten 1001 tot en met 10000. De toestand waarin de 1001ste klant het systeem aantreft zal voor iedere simulatierun anders zijn. Op deze manier worden n waarnemingen x1, …, xn verkregen welke onafhankelijk en representatief zijn voor het lange-termijn gedrag van het simuleerde systeem. Met (1) kan aldus een 95%-BI bepaald worden voor bijvoorbeeld de verwachte verblijftijd voor het voorbeeld. 2.
Batch Means methode
Bij de Batch Means methode wordt gebruikt gemaakt van één, maar daarentegen wel lange, simulatie run. De simulatie wordt opgedeeld in intervallen of subruns. Per interval wordt een gemiddelde bepaald. Voor bovenstaand voorbeeld zou x1 de gemiddelde verblijftijd van de eerste 10000 klanten zijn, x2 de gemiddelde verblijftijd van klant 10001 t/m 20000, etc. Op deze wijze worden n waarnemingen x1, …, xn verkregen die nagenoeg geen last hebben van het ‘initial bias’ effect, echter de waarnemingen zijn wel afhankelijk. Als bijvoorbeeld klant 10000 een heel druk systeem achterlaat, zullen de eerste klanten in het tweede interval ook lange verblijftijden hebben. Dit effect wordt uitgemiddeld als de intervallen of subruns voldoende lang zijn. Met (1) wordt vervolgens het 95%-BI bepaald. 3.
Random Replicatie methode
De Random Replicatie methode (RRM) werkt nagenoeg hetzelfde als de Replicatie methode, echter het stopmoment van elke run is iedere keer anders. In plaats van om de 10000 klanten een nieuwe gemiddelde waarneming te bepalen, stopt de simulatie na bijvoorbeeld 1 dag, 8 uur of 1 week. De n gemiddelde waarnemingen x1, …, xn zijn dus telkens gebaseerd op het aantal klanten dat in die periode aanwezig is geweest. Dit heeft als gevolg dat de gemiddelde waarnemingen niet meer identiek verdeeld zijn en dus ook (1) feitelijk niet meer gebruikt mag worden, althans niet voor de gemiddelde wachttijd per klant. Het volgende getallenvoorbeeld illustreert dit probleem van wissellende aantallen. Stel dat op de eerste dag 250 klanten binnenkomen en gemiddeld 6 minuten moeten wachten; op de tweede dag 150 klanten die gemiddeld 4 minuten wachten. Met de RRM zou dit een gemiddelde wachttijd opleveren van 5 minuten. De klanten hebben echter gemiddeld (240 x 6 + 150 x 4) / 400 = 5.25 minuten gewacht. M.a.w. drukke dagen worden even zwaar meegenomen als rustige dagen. Als de periode die in één simulatierun wordt gesimuleerd voldoende lang is (en de gemiddelden gebaseerd zijn op grote aantallen klanten), dan zal het effect van de wissellende aantallen minimaal zijn. In dat geval zou men met (1) een 95%-BI kunnen bepalen, maar statistisch gezien is dit niet 100% gerechtvaardigd.
