Frequentie en kans Tussen haakjes de namen van de programma's 1 2
3
4 5
6
7
(GRSIMUL1:) Gezond verstand? (KOPMUNT) ..............................................................................................................................1 Uitkomstenverzamelingen.........................................................................................................................................2 - Muntje gooien (KOPSTUK) ................................................................................................................................5 - Theorie en praktijk (GALTON, SITRINOM) ...................................................................................................5 - Combinaties ............................................................................................................................................................5 - Met twee dobbelstenen (DOB2SNEL) ...............................................................................................................9 - Drie dobbelstenen (DOB3PRAK, DOB3THEO) ................................................................................................9 - Meer munten (MUNTEN) .....................................................................................................................................10 - Pascal en Chevalier de Meré (PASCMERE ) ....................................................................................................11 - Het vaasmodel (VAASMODL ) .............................................................................................................................12 - Balletjes in bakjes leggen (BALINBAK ) .........................................................................................................13 - Series nullen en enen (SERIESIM ) ...................................................................................................................15 (GRSIMUL2 :) Het versterkte toeval.................................................................................................................................................16 - Van kop af winnen (AANKOP) .........................................................................................................................16 - Het vaasmodel van Pólya ( PROPOLYA ) ...........................................................................................................17 - Genetische drift (GENDRIFT ) ............................................................................................................................18 Wanneer verwacht je succes? .................................................................................................................................18 - Wanneer is de verzameling compleet? (COMPLEET ) .....................................................................................19 Beroemde problemen ...............................................................................................................................................21 - Het probleem van de twee kaarten (KAARTEN2 ) ............................................................................................21 - Het probleem van de oppas (OPPASSEN ) ........................................................................................................22 - Het probleem van de drie deuren (DEUREN3, DEURAUTO ) ......................................................................23 - Volledige inductie ..................................................................................................................................................24 - Pak twee getallen uit n. Hoe groot is de grootste van de twee? (MAX2UITN ) ..........................................25 - Eenentwintigen (TREFKANS ) .............................................................................................................................27 - Lootjes trekken voor Sinterklaas (SINT) ........................................................................................................29 - Schudden (SCHUDRAN , SCHUDIN) ....................................................................................................................30 - Matches (MATCHES) ............................................................................................................................................34 (GRSIMUL3 :) Continue kansvariabelen..........................................................................................................................................35 - Een duel (DUELPQ) ..............................................................................................................................................35 - Een meting van p, met darts (SIMULPI) .........................................................................................................35 - Ontwerp zelf een randomgenerator.....................................................................................................................26 - De afstand tot een cirkelrand (PINCIRKL ) .....................................................................................................29 - De afstand tot de rand van een vierkant (PIN4KANT , PINRHOEK ) .............................................................40 - De stomphoekige driehoek (AOBSTOMP ) .........................................................................................................41 - Een afspraakje (DATING) ...................................................................................................................................42 - Maak een driehoek van een stok (STOKIN3, STOKSNEL ) ............................................................................43 - De knekelman: het gebroken bot (STOKIN2, STOK2MED ) ...........................................................................45 - Onnodige ongerustheid (BAYES) ......................................................................................................................47 (GRSIMUL4 :) Extra opgaven over simulaties .............................................................................................................................50 (JAARSIM, RANDWALK , BIN4KANT , BINCIRKL , WORTELS2, STOWAND)
KANSEN en SIMULATIES 1
Gezond verstand?
De kansrekening kent veel verrassende uitkomsten. Het lijkt soms wel, of je te maken hebt met natuurkunde in plaats van wiskunde. Of erger: je gaat twijfelen aan "het" gezonde verstand. Dat zal je in het komende hoofdstuk vaak merken. De meest fundamentele controverse tussen theorie en praktijk is wellicht het experiment van de twee geldstukken. Werp met een euro. De kansen op kop of munt zijn ½, daarover kan geen twijfel bestaan. Werp tegelijk met twee identieke euro's en daar begint de eerste aarzeling. Wat is de kans op twee kop of twee munt? En wat is de kans op een kop en een munt? Het is heel goed mogelijk dat zelfs de meest onnozele beginneling het antwoord op deze vragen al weet, omdat hij of zij er al eerder over gehoord of gelezen had. Maar stel je eens voor dat je in de middeleeuwen leeft en de beschikking hebt over twee muntjes, een stuk papier en een ganzenveer en dat je over gezond verstand beschikt. Er zijn duidelijk twee antwoorden die overwogen moeten worden, samenhangend met het aantal mogelijkheden. En de hamvraag is of er drie dan wel vier mogelijkheden zijn. MODEL 1 K en K K en M M en M
MODEL 2 K en K K en M M en K M en M
Het eerste model is lijkt me het meest aantrekkelijk voor iemand met gezond verstand. De twee geldstukken immers zijn identiek, dus niet van elkaar te onderscheiden. De mogelijkheden K en M en M en K zijn dus ook niet van elkaar te onderscheiden, waarmee model 2 verworpen zou worden. Maar wacht eens. Als je de twee muntjes een verschillende kleur geeft of als je ze niet meer tegelijk, maar achter elkaar werpt, dan zijn ze in één klap wel van elkaar te onderscheiden! De rode van de groene of de eerstgeworpene van de tweede. Dat zou pleiten voor model 2. We gaan een keuze maken tussen de modellen en dat doen we met behulp van de randomgenerator A: randInt(0,1) . Een uitkomst 1 staat voor kop, het aantal keren kop wordt geteld in L2 : randInt(0,1)+randInt(0,1)üK 1+L2(K+1)üL2(K+1)
Onderaan het scherm verschijnt een teller. Omdat het werpen te snel gaat voor het menselijk oog, tellen we om de tien worpen: If
fPart(X/10)=0:Output(8,1,X).
Probeer KOPMUNT maar een paar keer uit met bijvoorbeeld 1000 worpen en je onzekerheid is voorgoed verdwenen: model 2 is het ware model!
KANSEN
-1-
ClrHome: Disp "HOEVEEL" Disp "SIMULATIES:" Input S ClrHome {0,1,2}üL1:3üdim(L2) For(X,1,S) randInt(0,1)+randInt(0,1)üK 1+L2(K+1)üL2(K+1) If fPart(X/10)=0:Output(8,1,X) End Output(1,1,"{K,K}:") Output(2,1,"{K,M}:") Output(3,1,"{M,M}:") Output(1,8,L2(1)) Output(2,8,L2(2)) Output(3,8,L2(3)) Output(5,1,"FRACTIE {K,M}:") Output(6,1,L2(2)/S) Output(8,7,"SIMULATIES")
De meest sceptische lezer vraagt zich zelfs af, of de random functie wel te vertrouwen is. Zijn de random waarden wel uniform (gelijkmatig) verdeeld tussen 0 en 1 en zijn opvolgende randomwaarden wel onafhankelijk van elkaar? Verderop zullen we hier dieper op ingaan. OPGAVEN 1.
De functie rand(999) levert 999 randomgetallen x tussen 0 en 1. De TI-83 kan er maar maximaal 999 produceren. Ga experimenteel na hoe groot x is, de gemiddelde waarde van deze randomgetallen. Bepaal het gemiddelde absolute verschil tussen twee randomgetallen:
| x - y | = abs (rand - rand ) ; is dat hetzelfde als | x - y | ? En is 2.
x´ y = x´ y
? Met andere woorden: zijn de uitkomsten onafhankelijk van elkaar?
INTERMEZZO
Een enquêteur ondervraagt 100 mensen over het bezit aan huisdieren. Uit zijn aantekeningen blijkt later dat 30 een hond hebben en 18 daarvan ook nog een ander huisdier 37 een poes hebben en 22 daarvan ook nog een ander huisdier 23 een konijn hebben en 10 daarvan ook nog een ander huisdier 4 mensen vertelden, van alle drie de soorten (hond, poes, konijn) er minstens een in huis te hebben. Er werden geen andere soorten gemeld. De enquêteur is helaas vergeten bij te houden hoe vaak het voorkwam dat er geen huisdier aanwezig was bij de ondervraagden. Help hem even!
2
Uitkomstenverzamelingen
Als alle mogelijke uitkomsten van een (toevals-)experiment bij elkaar worden gezet zonder op de volgorde daarvan te letten, krijg je de uitkomstenverzameling (U). Alle uitkomsten (elementen, elementaire gebeurtenissen) zijn verschillend, hoewel ze soms niet van elkaar te onderscheiden zijn. Dat is geen tegenspraak. Denk aan een bak met een witte en drie identieke rode knikkers. De rode knikkers zijn niet van elkaar te onderscheiden, maar voor de uitkomstenverzameling zijn het verschillende elementen. Het is gebruikelijk de elementen van een verzameling tussen accolades {} weer te geven. Je kunt U noteren als {W, R1, R2, R3} maar ook bijvoorbeeld U = {R2, W, R1, R3}, de volgorde is immers niet van belang.
KANSEN
-2-
Een ander voorbeeld. Bij het werpen met twee geldstukken, lettend op kop (K) of munt (M), is U = {(K,K), (M,M), (K,M), (M,K)}. De ronde haakjes geven aan dat daarbinnen de volgorde wél van belang is, zoals bij coördinaten: het punt (3,7) is een ander punt dan (7,3). Binnen accolades echter mogen de elementen door elkaar geschud worden. De ongeordendheid van een verzameling blijkt ook uit de definitie van gelijke verzamelingen. De verzamelingen A en B zijn gelijk aan elkaar (A=B) als elk element van A ook element van B is en omgekeerd elk element van B ook element van A is. Dit "in twee richtingen" toetsen van een eigenschap is kenmerkend voor bewijzen in de verzamelingenleer. Gebruikmakend van de symbolen Î voor "is een element van" of "ligt in" of "is lid van" en Þ voor "als … dan" staat er als definitie: A=B betekent: voor alle xÎA Þ xÎB en voor alle xÎB Þ xÎA. De "als dan en omgekeerd" symbolen kunnen samengetrokken worden tot Û hetgeen erg ondubbelzinnig maar typisch wiskundig wordt geformuleerd als "dan en slechts dan als": De verzamelingen A en B zijn gelijk dan en slechts dan als elk element van A element is van B en als elk element van B element van A is. Waaruit bijvoorbeeld volgt dat de verzamelingen {1,2,3} en {3,3,3,1,2} gelijk zijn. Enkele andere definities. A is een deelverzameling van B, notatie AÌB als voor alle xÎA Þ xÎB. Het complement van A, genoteerd als A', is de verzameling van alle elementen van U die niet tot A behoren. Kortweg geschreven als: A' = {xÎU | x Ï A}. De doorsnede van A en B is de verzameling van de elementen die in A en B liggen: AÇB = {xÎA | xÎB}; dit is hetzelfde als {xÎB | xÎA} want stel dat het niet zo was, dan zou er volgens de tweede definitie een element uit de doorsnede zijn dat niet in A ligt maar wel in A zou moeten liggen volgens de eerste definitie. Waarmee met een redenering uit het ongerijmde de gelijkwaardigheid is aangetoond. De vereniging AÈB is de verzameling van de elementen die in A of in B liggen met de aantekening dat het voegwoord of gelezen moet worden als en/of. De GR gebruikt hiervoor het woord or uit de logica (zie [2nd[TEST]). Dus AÈB = {x | xÎA or xÎB} De verschilverzameling A \ B of A – B is gedefinieerd als: {x | xÎA and xÏB)} = AÇB'. Het aantal verschillende elementen van A wordt genoteerd als #(A); de lege verzameling wordt genoteerd als Æ zodat #(Æ)=0.
KANSEN
-3-
Verzamelingen laten zich prettig illustreren met Venn-diagrammen. doorsnede AÇB and
leeg AÇB=Æ uitsluitend
vereniging AÈB or
complement A' not
verschil A\B
óf óf (A\B) È (B\A) xor
Dan nog een paar voorbeelden met getallen. Stel dat de uitkomstenverzameling U = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9} is en E = {0, 2, 4, 6, 8}; O = {1, 3, 5, 7, 9}; A = {0, 2} en B = {1, 2, 3, 4} dan geldt o.a.: EÌU, OÌU, AÌU, BÌU; AÌE; AÇB={2}, AÇE=A, AÇO=Æ, EÇO=Æ; AÈB={0, 1, 2, 3, 4}, AÈO={0, 1, 2, 3, 5, 7, 9}, EÈO=U; U\O = O' = E, U\E = E' = O, A\E =Æ, E\A = {4, 6, 8}; #(U)=10; #(E) + #(O) = # (EÈO). Deelverzamelingen A, B en C van een uitkomstenverzameling U hebben o.a. de volgende eigenschappen. PANORAMA
(1) (5) (7) (10) (12) (14) (16)
AÈÆ=A (2) AÇÆ=Æ AÈA=A (6) AÇA=A AÈA'=U (8) AÇA'=Æ AÈB = BÈA AÈ(BÈC) = (AÈB)ÈC AÈ(BÇC) = (AÈB)Ç(AÈC) (AÈB)' = A' Ç B'
(3) AÈU=U (9) (11) (13) (15) (17)
(4) AÇU=A
(A')' =A AÇB = BÇA AÇ(BÇC) = (AÇB)ÇC AÇ(BÈC) = (AÇB)È(AÇC) (AÇB)' = A' È B'
Deze stellingen zijn meestal gemakkelijk via een redenering te bewijzen. Daarbij moet vaak "heen" en "terug" een bewijs geleverd worden. Nemen we als voorbeeld stelling (15): AÇ(BÈC) = (AÇB)È(AÇC) Van links naar rechts. Een willekeurig element x van AÇ(BÈC) zal tot A en tot BÈC behoren, dus tot B of C. Dus x behoort tot A en B of tot A en C. Dus is x ook een element van (AÇB)È(AÇC).
KANSEN
-4-
Van rechts naar links. Stel y is een element van (AÇB)È(AÇC). Dan hoort y tot AÇB of tot AÇC. Als y element is van AÇB dan behoort y tot A en B dus ook tot BÈC. Als y een element is van AÇC dan kan met een soortgelijke redenering aangetoond worden dat y element is van A en BÈC. De twee verzamelingen links en rechts van het gelijkteken zijn dus identiek. Straks, bij het bewijs van de uitgebreide somregel, zullen we verwijzen naar enkele van deze stellingen.
Muntje gooien De eenvoudigste uitkomstenverzameling waar nog wat aan te beleven valt, is de verzameling {0, 1}. Je kunt daarbij denken aan de twee mogelijkheden bij het werpen met een geldstuk: bijvoorbeeld 0 = "munt" en 1 = "kop". {M, K} zou ook kunnen, maar we geven hier de voorkeur aan een uitdrukking met getallen. Met getallen immers kan je berekeningen uitvoeren, met de letters M en K gaat dat niet zo best. Bij een "ideale" munt is de verwachting dat munt en kop (0 en 1) gemiddeld ongeveer even vaak voor zullen komen. In het volgende programmaatje gaan we dat simuleren. De letter T telt de worpen (met 100 tegelijk), de letter F (van frequentie) houdt bij hoe vaak er een 1 werd geworpen. De breuk F/T heet de fractie (denk aan fPart). KOPSTUK ClrHome:ClrList L1 Disp "WERPEN MET EEN" Disp "GELDSTUK EN" Disp "TELLEN VAN DE" Disp "FRACTIE KOP" Output(8,14,"==>") Pause :ClrHome 0üT:0üF Output(1,1,"AANTAL: ") Output(2,1,"FRACTIE: ") Lbl 0 100+TüT randInt(0,1,100)üL1:sum(L1)üK F+KüF Output(2,10," ") Output(2,10,round(F/T,4)) Output(1,9,T) Goto 0
Als je dit programma uitprobeert, ben je een simulatie van een experiment aan het uitvoeren. Afrondend op drie decimalen ontdekte ik dat mijn TI-83 Silver Edition pas na een paar minuten min of meer tot rust kwam op de verwachte gemiddelde waarde 0.500; dat valt nog behoorlijk tegen. De vierde decimaal heeft veel langer nodig. Zie het schermplaatje: na ongeveer 7 minuten en 30 000 worpen zoekt de vierde decimaal nog steeds naar rust. Een supercomputer kan heel wat meer decimalen aan. Je kunt afschatten dat pas na ongeveer vier miljoen worpen (eindelijk!) de derde decimaal redelijk betrouwbaar is. [ Voor de gevorderden: normalcdf(2002000,E99,2000000,1000)=0.02275) ] Vier miljoen worpen: als er elke seconde een keer geworpen werd en men ging daarmee zeven weken lang door, dag en nacht, dan pas zou men de derde decimaal een beetje kunnen vertrouwen. De conclusie is dat we uit experimenten nooit met zekerheid de conclusie kunnen trekken dat de fractie uiteindelijk (op den duur) exact 0,5 is. Je kunt immers niet oneindig vaak gooien. Men moet dus zoeken naar een bewijs. Maar elk bewijs is gebaseerd op een aantal axioma's. Een axioma is een stelling die we niet (kunnen of willen) bewijzen, een soort principe dus. De Rus Kolmogorov (19031987) is er in geslaagd een bewijs te geven voor, wat we noemen, de wet van de grote aantallen.
KANSEN
-5-
Laten we in de volgende bladzijden Pierre-Simon Laplace (1749-1827) volgen, de Fransman die als eerste een model voor de waarschijnlijkheidsrekening opstelde. Na oneindig veel herhalingen zullen elementaire gebeurtenissen gemiddeld even vaak voorkomen. We veronderstellen dat de fractie dan overgaat in een kans. Een elementaire gebeurtenis (e) slaat op de elementen van de totaalverzameling. Uitgaande van een oneindige uitkomstenverzameling met n elementen: U = {e1, e2, e3, … , en} beschouwen we een deelverzameling A van U die we een gebeurtenis zullen noemen. Voor elke gebeurtenis AÌU definiëren we de kans als volgt:
P( A) =
# ( A) # (U )
Aangezien A een deelverzameling is van U en dus #(A) £ #(U) moet P tussen 0 en 1 liggen: 0 £ P(A) £ 1. Voor elkaar uitsluitende gebeurtenissen (verzamelingen met een lege doorsnede) geldt de somregel: P(AÈB) = P(A) + P(B) voor elk paar met AÇB=Æ We nemen deze eigenschap op als een soort axioma en bewijzen hem niet. De uitgebreide somregel echter, die ook geldig is voor niet gescheiden verzamelingen: P(AÈB) = P(A) + P(B) - P(AÇB) moet bewezen worden. En dat doen we met behulp van enkele hiervoor gegeven stellingen uit de verzamelingenleer. Het bewijs verloopt als volgt. (*)
AÈB = AÈ(B\A), immers: AÈB = (AÈB)ÇU (stelling 3) = (AÈB)Ç(AÈA') (stelling 7) = AÈ(BÇA') (stelling 14) = AÈ(B\A)
(**)
(B\A) È (BÇA) = B, want: (BÇA') È (BÇA) = B Ç (AÈA') (stelling 15) = B Ç U (stelling 7) = B (stelling 3)
Wat de kansen betreft geldt dus volgens het axioma en volgens (**): P(B\A) + P(AÇB) = P(B), want de doorsnede is leeg en AÇB = BÇA. Waarna uit (*) volgt dat P(AÈB) = P(AÈ(B\A)) = P(A) + P(B\A) = P(A) + P(B) - P(AÇB). Zonder in te gaan op de finesses (dat valt buiten het kader van dit boek) vermelden we de uiterst belangrijke productregel. Werp je bijvoorbeeld twee keer met een dobbelsteen dan mag aangenomen worden dat de uitkomst van de eerste worp niet van invloed is op de tweede worp. Omdat bij elk van de zes mogelijkheden van de eerste worp weer zes mogelijkheden voor de tweede worp horen, zijn er in totaal 6´6=36 mogelijkheden. De kans op een zes bij elke worp is 1/6. De kans om twee achter 1 1 1 elkaar een zes te gooien is dus = ´ ; waarmee de productregel geïllustreerd is: 36 6 6 Voor de onafhankelijke gebeurtenissen A en B geldt: P(A en B) = P(A)´P(B). Laten we ook dit principe als een gegeven accepteren. OPGAVE 3.
Creëer 400 randomgetallen in L1, 400 nieuwe randomgetallen in L2 en doe L1*L2->L3. Hoe kun je nu via de gemiddelde waarden de (on)afhankelijkheid van al deze random trekkingen experimenteel onderzoeken?
KANSEN
-6-
Combinaties Op hoeveel manieren kun je van (0,0) naar (2,4) bewegen in gehele stapjes naar rechts (positieve X-richting) of omhoog (positieve Y-richting)? Getekend is de route (0,0) ® (0,1) ® (0,2) ® (1,2) ® (1,3) ® (1,4) ® (2,4). De letters X en Y gebruikend kan dit sneller genoteerd worden als YYXYYX en nu is de vraag, hoeveel rangschikkingen er mogelijk zijn van de letters XX YYYY. Met genummerde letters staat er: X1X2 Y1Y2Y3Y4. Op de eerste plek kan gekozen worden uit 6 tekens, bij elk van die keuzes zijn er 5 mogelijkheden voor de volgende plek en zo door. In totaal dus 6.5.4.3.2.1 = 6! rangschikkingen. Echter zijn veel van die mogelijkheden hetzelfde omdat de letters niet genummerd zijn. Zo is Y1X1X2Y4Y2Y3 niet te onderscheiden van Y3X2X1Y4Y2Y1. De 2 letters X zijn op 2! manieren te ordenen, de 4 letters Y op 4! manieren. Het totaal aantal routes van (0,0) naar (2,4) is dus 6! 6.5.4.3.2.1 = = 15 ; 2! 4! 2.1.4.3.2.1
dit noemt men de combinaties van 2 uit 6. Op de TI-83 worden de combinaties in dit geval berekend met MATH
6 nCr 2. Stel dat de kans op een stapje naar rechts p is en de kans op een stapje omhoog q (dus p+q=1) dan is de waarschijnlijkheid om vanuit (0,0) in (2,4) terecht te komen: P (2,4) =
6! 2 4 . p .q ; 2! 4!
zo ontstaat de binomiale kansverdeling. In het geval dat p = q = 0,5 is de uitkomst P ( 2,4) = 15.0,56 =
15 64
We gaan nu over op een toevalswandeling. Op de lijn x + y = 6 liggen de roosterpunten waarvoor in totaal 6 keer een aselecte (random) keuze tussen X of Y moet worden gemaakt. De aantallen routes staan hier aangegeven, de bijbehorende kansen zijn: P(0,6)=P(6,0)= 641 ; P(1,5)=P(5,1)= 646 ; 20 P(2,4)=P(4,2)= 15 64 ; P(3,3)= 64
Dit binomiale kansproces is op een bekende manier zichtbaar te maken met het spijkerbord van Galton. Maak spijkers op de plaats van de roosterpunten, draai het bord zo, dat het spijkertje (0,0) bovenaan komt en laat een groot aantal knikkers van bovenaf naar beneden stuiteren. Als alles mooi symmetrisch en horizontaal is opgesteld is de kans op linksaf even groot als de kans op rechtsaf en ontstaat een ongeveer symmetrisch histogram patroon. In het vorige plaatje kun je de dikke stippen op de lijn beschouwen als de verzamelbakjes van het spijkerbord. Kijk even naar het (helaas vierkante) knikkertje dat vlak boven de tweede spijker zweeft en moet kiezen tussen linksaf of rechts af. Een knikker die 6 keer linksaf geslagen is (linksaf, gezien door de ogen van de lezer, rechtsaf vanuit het gezichtspunt van het knikkertje) komt in opvangbakje nummer 0 terecht. De broncode van het simulatieprogramma GALTON dat de binomiale verdeling illustreert laten we achterwege. Er is een vertragingsfactor U ingebouwd (anders kun je de vallende knikkers nauwelijks met het blote oog volgen). De kernregels van het programma zijn gemarkeerd.
KANSEN
-7-
OPGAVE 4.
Vergelijk de hierboven gevonden aantallen met de theoretisch verwachte aantallen. Gebruik binompdf.
De driedimensionale uitbreiding van het binomiale kansexperiment heet trinomiaal. Het bijbehorende simulatieprogramma heet SITRINOM. Bekijken we als voorbeeld de wandeling van (0,0,0) naar (4,5,3) via 2 stappen omhoog, 4 stappen naar voren, 3 naar rechts, 1 omhoog en 2 naar rechts. Symbolisch: 12 stappen in de volgorde ZZXXXXYYYZYY. Hoeveel verschillende wegen leiden naar het punt (4,5,3)? Kies er één speciale volgorde uit bijvoorbeeld: XXXX YYYYY ZZ. Net als in het tweedimensionale geval realiseren we ons dat het totaal aantal rangschikkingen gelijk is aan 12!, dat de X'en op 4! manieren, de Y'en op 5! manieren en de twee Z'en op 2! manieren gerangschikt kunnen worden zonder dat dat tot een ander resultaat zou leiden. De conclusie moet zijn dat het aantal routes gelijk is 12! aan . 4! 5! 3! Zijn de kansen op een stapje in de x-, y- en z-richting respectievelijk p, q en r (met p+q+r=1) dan kan de kans dat een knikkertje in een denkbeeldig driedimensionaal spijkerbord in het vakje (4,5,3) terechtkomt geschreven worden als: 12! . p 4 .q 5 .r 3 . 4! 5! 3! Uitgaand van gelijke kansen in de drie richtingen wordt dit: P (4,5,3) =
12! 1 4 1 5 1 3 . . . » 0,05216 4! 5! 3! 3 3 3
Voor de controle van deze theorie via een simulatie moeten we een kleiner aantal stappen nemen, om te voorkomen dat ons eenvoudige rekenapparaat urenlang moet blijven rekenen. We gaan uit van slechts vier stappen en gelijke kansen in de drie richtingen. Ons virtuele driedimensionale spijkerbord is dan de gelijkzijdige driehoek die hier getekend is in het vlak met vergelijking x+y+z=4. De dikke stippen stellen weer de opvangbakjes voor. In het binnengebied van de driehoek liggen de bakjes (2,1,1), (1,2,1) en (1,1,2). Het aantal wegen 4! naar die punten is = 12 en de kans daarop is 1!1!2! 12 4 12. 1 3 = » 0,148 . 81 4
Om gehele verwachtingswaarden te krijgen is een veelvoud van 3 (dus van 81) handig. Gekozen is voor 810 simulaties. De verwachting is dan dat de drie aangegeven punten {1,1,2} = ((1,1,2), (1,2,1), (2,1,1)} ongeveer 3´120 keer een knikker opvangen. De punten aan de buitenrand van de driehoek zullen naar verwachting ongeveer 3´10, 6´40 en 3´60 keer geraakt worden. Het experiment levert dit keer een verbluffend resultaat (rechts), haast te mooi om waar te zijn. Als je zo'n uitslag ziet, wees dan op je hoede! Inderdaad, tegen mijn gewoonte heb ik dit keer gesjoemeld (net zolang uitgeprobeerd tot het resultaat me beviel).
KANSEN
-8-
Met twee dobbelstenen (kijkend naar de som van de ogen) Deze simulatie is al eerder aangestipt. We geven nu de kern van een programma dat sneller werkt. Het nadeel daarvan is dat je tussentijds de uitkomsten van de worpen niet kunt zien. Er wordt namelijk 500 keer tegelijk geworpen met de twee dobbelstenen. Eén ronde van 500 worpen kost ongeveer 16 seconden op de TI-83 SE (de gewone TI-83 werkt ongeveer twee keer zo traag). Het voordeel van deze methode is dat je binnen een redelijk tijdsbestek een erg fraai resultaat krijgt. De theorie. Op hoeveel manieren krijg je een som acht? De "gunstige" gebeurtenissen met X+Y=8 ( (6,2) t/m (2,6) ) kunnen we aangeven in een roosterpatroon van 6 bij 6. Zo wordt, bijvoorbeeld, de kansberekening P(X+Y=8) geïllustreerd en het antwoord 5/36 onmiddellijk begrepen. We verwachten dus ongeveer 3000×5/36 = 417 keer som acht na 3000 worpen. DOB2SNEL seq(X,X,1,12)üL1 For(J,1,N) randInt(1,6,500)+randInt(1,6,500)üL3 For(I,1,12) L2(I)+sum(L3=I)üL2(I) End End
In de vijfde regel wordt er, iedere keer als L3=I waar is (dus de waarde ervan 1 is) 1 bij sum(L3) opgeteld. Aldus sommeer je de gunstige uitkomsten in L2(I).
Drie dobbelstenen Het toevoegen van een enkele dobbelsteen compliceert de theoretische zaak aanmerkelijk. Kijken we bijvoorbeeld naar het geval dat het ogentotaal van de drie dobbelstenen negen is: x+y+z=9 {1,2,6} heeft {1,3,5} heeft {1,4,4} heeft {2,2,5} heeft {2,3,4} heeft {3,3,3} heeft Totaal:
KANSEN
3!= 3!= 3!=
6 rangschikkingen 6 rangschikkingen 3 rangschikkingen 3 rangschikkingen 6 rangschikkingen 1 rangschikking ---- + 25 rangschikkingen
-9-
Het totaal aantal mogelijkheden is 6´6´6=216 dus P(x+y+z=9) = 25/216 = 0.11574 Het is erg moeilijk om daar een model voor te vinden. Tot en met een som s=8 verloopt het aantal mogelijkheden volgens de rij ½(s-1)(s-2) maar daarna klopt die formule niet meer. Echter ... Wie niet slim is moet sterk zijn. We laten de rekenmachine het werk voor ons doen. Hij krijgt de opdracht domweg alle mogelijkheden te onderzoeken in drie "geneste" for-lussen. Zie het programmadeel van DOB3THEO. S is de som van de ogen, T is de teller. In L2 komen de verwachte frequenties, voor de kansverdeling moet je L2 nog door 216 delen. DOB3THEO 0üT For(S,3,18) For(X,1,6) For(Y,1,6) For(Z,1,6) If X+Y+Z=S:T+1üT End:End:End Output(1,5,S) TüL2(S-2):0üT End
Het programma DOB3PRAK is gemakkelijk te construeren, na DOB2SNEL. Ik geef alleen het fundament en de schermplaatjes na 50´216 simulaties. For(J,1,N) randInt(1,6,216)+randInt(1,6,216)+randInt(1,6,216)üL3 For(I,3,18) L2(I-2)+sum(L3=I)üL2(I-2) End End
Meer munten Het programma MUNTEN simuleert het werpen met maximaal vijf geldstukken. Je kunt de teller "aantal kop" zien lopen tijdens de simulatie.
KANSEN
- 10 -
Pascal en Chevalier de Meré Het meningsverschil tussen Pascal en Chevalier de Meré. Ridder Chevalier de Méré (eind 17de eeuw) wedde met succes dat hij in 4 worpen met een dobbelsteen minstens één zes zou gooien (spel A). De kans daarop is namelijk ietsje groter dan 0,5 (ongeveer 0,5177) terwijl zijn simpele verwachting 4. 1 6 = 2 3 was. Ook probeerde hij spel B uit: na 24 worpen met twee dobbelstenen hoopte hij minstens één dubbelzes te gooien met dezelfde kans (volgens dezelfde simpele redenering: 24. 1 36 = 2 3 ). Het verhaal gaat dat tot de Méré's teleurstelling en woede (na uitgebreid experimenteren) bleek dat de kans op succes in spel B iets kleiner dan 0,5 was. De wiskundige Pascal toonde dit aan en ontwikkelde zo een van de eerste kansmodellen. De kansen op minstens één zes kunnen we bepalen met de complement regel. In TI-83 taal: 1 - binompdf(4,1/6,0) » .5177 voor A en 1 - binompdf(24,1/36,0) » .4914 voor B
Ook met de kans op precies één (dubbel)zes kom je niet op ongeveer ½ en al helemaal niet in de buurt van 2 3 uit: binompdf(4,1/6,1) » .3858: binompdf(24,1/36,1) » .3488 Als je het eindeloze werpen van De Meré simuleert, begint de conclusie zich op te dringen dat het hoogst onwaarschijnlijk is dat hij voldoende tijd en uithoudingsvermogen had om het kleine verschil tussen kansA en kansB experimenteel aan te tonen. Iemand heeft de geschiedenis dus een beetje aangedikt. Een schatting. Hoe vaak moet je (bij kansB = 0,4914) 24 keer met twee dobbelstenen gooien, om met (zeg) 95% zekerheid te kunnen constateren dat het gemiddelde m kleiner is dan 0,5? Gebruikmakend van m = np en s = Ö(np(1-p)) » Ö(n.½.½) = ½Ön en uitgaand van z=invNorm(0.05)=-1.645 komt er volgens de regels van het toetsen: (0,5n-0,4914n) = 1,645.0,5Ön waarvan de oplossing is n>9091. [Of, geheel in TI stijl via MATH 0: 0=normalcdf(-E99,.4914N,.5N,.5Ö(N))-.05 ] Alleen al om met redelijke zekerheid te constateren dat kansB<0,5 is moet je dus minstens 9000 keer gooien; en dan is er nog dat andere experiment (kansA=0,5177)! Stel dat De Meré 10 keer per minuut wierp met een paar dobbelstenen (hij was ook tijd kwijt met het noteren, met een ganzenveer, ballpoints zijn pas na 1960 in zwang gekomen) en dat 24 keer, dan moet hij voor 9000 experimenten in totaal minstens 40 uur bezig geweest zijn. Waar verslaving al niet toe kan leiden... OPGAVE 5.
Reken uit hoeveel keer er minstens geworpen moet worden met een geldstuk, om met 95% waarschijnlijkheid te kunnen besluiten tot p = 0,50 op twee decimalen nauwkeurig. Ga dus uit van 0.495 < p < 0.505.
PASCMERE For(I,1,100) randBin(4,1/6)üP If P>0:A+1üA randBin(24,1/36)üQ If Q>0:B+1üB Output(1,9,I) Output(4,4,A) Output(5,4,B) End
KANSEN
- 11 -
Het vaasmodel Eerst de theorie achter het experiment in VAASMODL. Gegeven: in een vaas liggen z zwarte en w witte knikkers (totaal z+w=n). Hieruit wordt m keer aselect een knikker gepakt: daarvan zijn x zwart. æ mö z w MET terugleg is de kans daarop çç ÷÷. p x .q m - x met p = en q = z+w z+w è xø Dit is het binomiale model binompdf(m,p).
æzöæ w ö çç ÷÷.çç ÷ x ø è m - x ÷ø è ZONDER terugleg is de kans ; het hypergeometrische model. æ z + wö çç ÷÷ è m ø Z zwarte en W witte knikkers (samen N knikkers) zitten in de vaas; er worden P knikkers uit gepakt. X is het aantal getrokken zwarte knikkers (het aantal succesen); X loopt van 0 tot P. Er zijn S simulaties.
Hier volgt de simulatie waarbij getrokken wordt MET teruglegging. Het aantal getrokken zwarte knikkers komt in L2; de theoretische kansverdeling komt in L3. For(J,1,S) 0üX For(I,1,P) randInt(1,N)üR If randInt(1,N)÷Z:X+1üX End L2(X+1)+1üL2(X+1) Output(8,1,J) End round(S*binompdf(P,Z/N),0)üL3
KANSEN
- 12 -
En hier komt de trekking ZONDER teruglegging. XüN:YüM For(J,1,S) For(I,1,P) If randInt(1,N+M)÷N
"ZWARTE GEPAKT ZONDER TERUGLEG"
Then N-1üN:T+1üT Else M-1üM End End L2(T+1)+1üL2(T+1) XüN:YüM:0üT Output(8,1,J) End For(K,0,P) (X nCr K)((Y nCr (P-K))/((X+Y) nCr P)üQ round(SQ,0)üL3(K+1) End
"THEORIE"
OPGAVEN 6.
Trek n knikkers tegelijk uit een vaas met 1 zwarte en 99 witte knikkers. Bereken de kans dat de zwarte knikker er bij zit. Simuleer dit met n=37.
7.
In één van 10 dozen zit een geldprijs. Tien deelnemers moeten om de beurt een doos openen, tot iemand de prijs gewonnen heeft. Wat is de gunstigste positie: als eerste deelnemer een doos openen of als laatste? Of maakt het niet uit? Simuleer de keuze dat je als laatste een doos wilt openen met het programma VAASMODL. Wat is (na bijvoorbeeld 100 simulaties) de kans dat je 9 voorgangers allen gefaald hebben zodat de prijs in het laatste, door jou nog te openen doosje zit?
Balletjes in bakjes leggen Bekeken wordt de kansverdeling van M balletjes die in N bakjes worden geplaatst. Het spiegelbeeld van trekken met terugleg is het plaatsen zonder terugnemen. Het wekt dan ook geen verwondering dat we in het laatste geval te maken hebben met een binomiale verdeling met de inverse kans P = 1/N, als N het aantal bakjes is waarin de balletjes geplaatst worden. Neem, om dit verder te verduidelijken, bijvoorbeeld vier balletjes (M=4) die aselect in twee bakjes (N=2) geplaatst worden. X is het aantal balletjes in een willekeurig bakje. Hoe groot is de kans dat er respectievelijk 0, 1, 2, 3, 4 balletjes in terechtkomen? P(X=0) is de kans dat er geen balletje in komt. 4 Die kans is ½.½.½.½ = (½) = 1. 116 P(X=1) vinden we door te bedenken dat één van de vier balletjes uitverkoren wordt en de andere drie niet. Er zijn vier mogelijkheden om uitverkoren te worden (want er zijn vier balletjes). Dus P(X=1) = 4. 1 2 .( 1 2 ) 3 = 4. 116 Vervolgens. Hoeveel combinaties van twee uitverkoren balletjes uit vier zijn er? Het woord verraadt 4! de uitkomst. Er zijn = 6 combinaties. 2!.2! Dus P(X=2) = 6. 1 2 2 . 1 2 2 = 6. 116 ; Verder nog P(X=3) = 4. 116 en P(X=4) = 1. 116
KANSEN
- 13 -
1
234
12
34
2
134
13
24
3
124
14
23
4
123
23
14
24
13
34
12
1234
Herken de binomiale verdeling! Als het gaat om M balletjes in N vakjes dan is de succeskans p = 1/N (want alle bakjes hebben een even grote kans, ze worden immers aselect gevuld) en er moeten X combinaties gemaakt worden (X = 0, 1, 2, …, M). De verwachtingswaarde van X wordt dus berekend met binompdf(M,1/N,X). BALINBAK 0üP ClrAllLists seq(X,X,1,N)üL1 For(K,1,S) ClrList L2:Nüdim(L2) For(X,1,M) randInt(1,N)üJ L2(J)+1üL2(J) End max(L2)üQ If Q>P:QüP P+1üdim(L3) P+1üdim(L„) seq(X,X,0,P)üL3 For(L,0,P) L„(L+1)+sum(L2=L)üL„(L+1) End Output(6,12,K) End
"S SIMULATIES" "BAL NR.X IN BAKJE J OPGETELD"
"IN L4 DE BAKJES MET L BALLEN" "K = SIMULATIENR."
Pause:ClrHome Disp L3 ,L 4 Pause
De laatste regels zijn toegevoegd om tussentijds te laten zien hoe vaak er 0, 1, 2, 3 of 4 balletjes in een bakje komen. Op het schermplaatjes hierboven is dat te zien na 100 simulaties. Hierna wordt aangegeven in L2 welke fractie van de ballen in een willekeurig bakje terechtgekomen is. Hoe fraai de fracties in L2 (altijd volgens het toeval) kloppen met de theoretische kansverdeling in L3 zie je bij de rest van het programma:
KANSEN
- 14 -
L„/(SN)üL„ For(X,0,P) binompdf(M,1/N,X)üL…(X+1) End ClrList L1,L2 L3üL1 :L„üL2:L…üL3 ClrList L„,L…:ClrHome
"FRACTIE" "THEORIE"
Als er 50 balletjes over 100 vakken verdeeld worden ontstaat (omdat n groot is en p klein) bij benadering een Poissonverdeling met p=0.01 en n=50 dus gemiddelde m=0.5. In formule: P ( x = k ) =
m k -m ×e k!
Series nullen en enen Als een rij enen en nullen (denk aan worpen met een geldstuk) gesimuleerd wordt, valt op dat er nogal wat lange rijen gelijke uitkomsten (series) bij zijn. Het volgende programma telt de series met gelijke lengte. De kans op een L+1 serie met lengte L is ½ . Dat is eenvoudig in te zien. Een serie met L=5 verschijnt als MKKKKKM of KMMMMMK. Noem de i-de uitkomst u i, dan is bij een serie met lengte L: P(u i-1 ¹ u i) = ½ ; P(u i = u i+1) = ½ ; P(u i+1 = u i+2) = ½ … P(u i+L-1 ¹ u i+L) = ½. L+1 Uitgaande van de onafhankelijkheid van alle worpen krijg je zo het product ½ . Hierbij is geen L rekening gehouden met de randen van de hoofdrij. De eerste serie namelijk heeft kans ½ en het laatste stuk van de rij kan geen serie meer bevatten. Vandaar dat we de series halen uit een oneindige rij. Aangezien de lijsten veel geheugenruimte opslokken, beperken we ons tot 512 simulaties die tijdelijk in L3 worden opgeslagen: randInt(0,1,512)üL3. Begin met serielengte L=1. De voorafgaande redenering volgend, wordt gezocht naar opvolgend gelijke uitkomsten: If L3(X)=L3(X+1)
Als dat het geval is, wordt de lengte van de serie (L) met 1 verhoogd: L+1üL Zoniet, dan wordt L weer op 1 gezet. Aan het eind wordt L3 schoongeveegd, om plaats te maken voor L+1 de theoretische verwachting 512 . ½ : {128,64,32,16,8,4,2,1,.5,.25,.125}üL3
KANSEN
- 15 -
SERIESIM 1üL seq(X,X,1,15)üL1 Output(7,1,"WORPNR.:") For(X,1,511) Output(8,1,X) If L3(X)=L3(X+1) Then L+1üL Else L2(L)+1üL2(L) 1üL End End L2(L)+1üL2(L)
Omdat het turven nogal snel gaat, is het verleidelijk om meerdere runs van 512 simulaties uit te voeren en pas daarna naar het totaal gemiddelde te kijken. Volgens de wet van de grote getallen lijken de frequenties steeds meer op de door de theorie verwachte waarden. Je kunt voor 100 runs van 512 worpen beginnen met: ClrAllLists For(R,1,100) 512üdim(L3 )
…. End
waarna de frequenties gecumuleerd (gestapeld) in L2 komen. De theoretische waarden moeten dan met 100 vermenigvuldigd worden. In de eerste sessie van 512 worpen kwam een recordlengte voor van 14 enen en dat record werd niet meer verbeterd in de negenennegentig sessies daarna.
3
Het versterkte toeval
Van kop af winnen Met toevalsprocessen gaat het anders dan bij de eindsprint van een wieleretappe. Wie in de statistiek aan kop ligt, blijft aan kop (in het algemeen). Stel we werpen twintig keer met een muntje. Als het gemiddelde begint met de waarde 1 (eerste worp) dan zal in het algemeen het voortschrijdend gemiddelde van bovenaf de ½ naderen. Hoe komt dat? Veronderstel dat je eerste worp een 1 was, gemiddelde dus 1. Na de tweede worp is de verwachting van het gemiddelde gedaald tot het gemiddelde van 1 en ½ is 0,75; na de derde worp tot 2/3 » 0,67. Bekijken we het gemiddelde na 4 worpen, dan is de verwachting dat het gemiddelde 0,625 zal zijn. Kijk maar: eerste tweede derde vierde worp 1 1 1 1 1 1 1 0 1 1 0 1 1 1 0 0 1 0 1 1 1 0 1 0 1 0 0 1 1 0 0 0
gemiddelde 1 0,75 0,75 0,5 0,75 0,5 0,5 0,25
Totaal gemiddelde 0,625 (optellen en door acht delen).
KANSEN
- 16 -
AANKOP For(X,1,S) Output(3,12,X) randInt(0,1,20)üL1 cumSum(L1=0)üL2 cumSum(L1=1)üL3 sum(L2ùL3)üI L„(I+1)+1üL„(I+1) End
Het vaasmodel van Pólya Het vaasmodel (urn model) van Pólya beschrijft het volgende. In een vaas zit één rode en één witte knikker. Pak er een uit, leg die weer terug en voeg een knikker van dezelfde kleur toe. Herhaal dit. Wat gebeurt er met de verhouding rood/wit in de vaas na een groot aantal trekkingen? Je zou drie eindsituaties kunnen verwachten: (1) ten slotte zijn er vrijwel 100% rode (of witte) knikkers in de vaas (2) de evenwichtssituatie zit op één waarde tussen de 0 en 100% (bijv 50%-50%) (3) de situatie blijft instabiel Verrassend genoeg voorspelt geen van deze drie opties de werkelijkheid! De vierde, juiste optie is: (4) Er komt uiteindelijk een evenwicht, maar je kunt niet voorspellen in welke verhouding rood/totaal het evenwicht bereikt wordt. Dit Polya proces komt heel wat vaker voor in de praktijk dan je zou denken. Bij koersschommelingen op de beurs bijvoorbeeld, of bij de natuurlijke selectie (genetische drift) blijkt dat kleine toevalstreffers in het begin van grote invloed zijn op het eindresultaat door dit principe van de positieve, versterkte terugkoppeling. Een bekend voorbeeld is dat van de eerste video-recorders. Aanvankelijk waren er gelijktijdig twee gelijkwaardige en even dure types (VHS en Betamax). Om onverklaarbare redenen (beter: door onverklaarbare oorzaken) nam het type VHS in het begin een kleine voorsprong, waarna Betamax uit de markt verdween. Niet door een kwaliteitsverschil of prijsverschil, maar door het versterkte toeval.
PROPOLYA Lbl 0 AüR:BüW:R+WüN For(J,1,100) For(I,1,P) R/NüS If rand<S Then R+2üR Else W+2üW End N+2üN End Pt-On(J,S) End Goto 0
KANSEN
- 17 -
Genetische drift Als een kleine populatie zich voortplant in een geïsoleerde situatie (op een eiland bijvoorbeeld), zullen kleine genetische veranderingen (mutaties) zich versterkt en versneld ontwikkelen. Onder die voorwaarden is de vorming van nieuwe soorten waarschijnlijker dan in gewone omstandigheden. Dat heet genetische drift. We gaan in het programma GENDRIFT uit van een startwaarde p = 0,5 en veranderen die geleidelijk aan door bij iedere stap (X+1->X stelt een nieuwe generatie voor) een verandering in p aan te brengen (volgens binomiaal toeval) die omgekeerd evenredig is met de grootte van de populatie (N): randBin(N,P)/N. Zie hoe na 100 generaties bij een kleine startpopulatie N=3 de eindtoestanden ver uiteen liggen (de nieuwe soorten); en zie hoe bij de relatief grote populatie N=20 de eindtoestand zich stabiliseert rond één constante p-waarde » 0,5. Vergelijk deze onverwachte en verrassende plaatjes met die van de twee vorige programma's AANKOP en PROPOLYA. De kern van dit programma is simpel: GENDRIFT Lbl A randBin(N,P)üT X+1üX:S+T/NüS:S/XüY line(A,B,X,Y) YüP:XüA:YüB If X<100:Goto A
4
Wanneer verwacht je succes?
In een vaas liggen 10 knikkers, 8 rode en 2 witte: {R,R,R,R,R,R,R,R,W,W}. Hieruit wordt aselect een aantal keren een knikker gepakt met teruglegging. De kans op een witte knikker is p, de kans op een rode knikker is q = 1-p. In dit voorbeeld dus p = 0,2 en q = 0,8. De verwachtingswaarde van een gebeurtenis is al eerder ter sprake gekomen. Daarvoor moet je alle voorkomende waarden met hun kans vermenigvuldigen en optellen: S(xi.pi). Na hoeveel trekkingen (gemiddeld) verwachten we voor het eerst een witte knikker? Het aantal pogingen (trekkingen) noemen we xi = 1, 2, 3 … De verwachtingswaarde van het gemiddeld aantal pogingen X is als volgt opgebouwd. Als je direct succes hebt (na de eerste poging al een witte getrokken hebt) is xi.pi gelijk aan 1´0,2; voor het daarop volgende succes zijn er X+1 pogingen nodig (dus: de eerste keer nog geen succes, daarna verloopt het volgens de verwachting X pogingen later) en komt er 0,8.(X+1) bij. In formule: X = S(xi.pi) = 1.0,2 + 0,8.(X+1) waaruit X opgelost kan worden: X = 0,2 + 0,8X + 0,8 dus 0,2X = 1 en X =
1 = 5. 0,2
In het algemeen (met een succeskans p in plaats van 0,2 en q = 1-p) krijg je X = 1.p + q.(X+1) uitgewerkt tot X = p + qX + q en (1-q)X = p+q. 1 Omdat p+q=1 is kun je hiervoor schrijven: p.X = 1 en X = p
KANSEN
- 18 -
De verwachtingswaarde van het gemiddeld aantal pogingen is dus het omgekeerde van de succeskans. Dat is natuurlijk nog zo gek niet. Als je met een dobbelsteen gooit (kans op een zes (succes) is 1/6), dan is het te verwachten aantal worpen voordat je 1 een zes gooit = 6 . Dat controleren we met de TI-83 via: 1
6
sum(randInt(1,6,999)=6)/999. Het gemiddeld aantal pogingen is dus
ongeveer 0,16 à 0,17: goed in overeenstemming met 1/6 » 0,167. OPGAVE 8.
Terug naar de vaas met tien knikkers. Noemen we het trekken van een rode knikker succes (p = 0,8) dan geldt evenzeer dat het verwachte aantal trekkingen 1/p is, dus verwachten wij gemiddeld 1/0,8 = 1,25 pogingen voordat een rode knikker getrokken wordt. Controleer dit een aantal keren met rand(999)<.8. Met rand<.8 worden toevalsgetallen tussen 0 en 0,8 geproduceerd om het trekken van een rode knikker te simuleren.
Wanneer is de verzameling compleet? Een dobbelsteen wordt geworpen. Na hoeveel pogingen (gemiddeld) zijn alle zes de mogelijkheden aan de beurt geweest? Met andere woorden: wanneer is de verzameling {1, 2, 3, 4, 5, 6} compleet? Gebruikmakend van het voorafgaande kunnen we het volgende stellen. Noem een geschikte worp een worp die nog niet aan de beurt geweest is. Noem de wachttijd van een worp het gemiddeld aantal nodige worpen voordat hij geschikt is. Worp nummer 1 is altijd nodig; worp nummer 2 is geschikt met een kans p = 5/6 dus de wachttijd is 1/p = 6/5; worp nummer 3 is geschikt met een kans p = 4/6 dus de wachttijd is 1/p = 6/4; enzovoorts. Het gemiddeld aantal verwachte worpen voor een complete verzameling is dus: 1+
6 6 6 6 6 6 6 6 6 6 6 1 1 1 1 1 1 + + + + = + + + + + = 6.( + + + + + ) 5 4 3 2 1 6 5 4 3 2 1 1 2 3 4 5 6
met de GR berekend via 6*sum(seq(1/X,X,1,6))=14.7 In theorie geldt dus voor het verwachte aantal worpen om een verzameling van n elementen compleet te krijgen: n.(1 + 1 2 + 1 3 + 1 4 + ... + 1 n ) . De controle op de theorie wordt uitgevoerd in het programma COMPLEET, in de lus: Repeat I=N (herhaal tot I gelijk is aan N) waarbij telkens een element A random uit de uitkomstenverzameling gekozen wordt. sum(A=L1)=0 onderzoekt of A al in L1 stond. Zo nee, dan wordt A alsnog in L1 geplaatst. X telt het aantal benodigde pogingen bij elke simulatie.
KANSEN
- 19 -
OPGAVE 9.
Het gemiddeld aantal trekkingen om een serie compleet te krijgen valt nogal mee. a) Na hoeveel trekkingen (gemiddeld) heb je een serie van 100 compleet? b) Na gemiddeld 2000 trekkingen is een bepaalde serie van n compleet. Bereken n met de MATH 0:Solver en controleer je oplossing met COMPLEET .
COMPLEET Lbl A ClrList L1 0üI:0üX "X=AANTAL POGINGEN" 0üL1(1) "L1 WORDT OPGEBOUWD" Repeat I=N "HERHAAL TOT I=N" randInt(1,N)üA "A=ELEMENT VAN UITKOMSTENVERZ" If sum(A=L1)=0 "ALS A NOG NIET IN L1 ZIT" Then I+1üI:AüL1(I) "DAN WORDT A IN L1 GEPLAATST" End X+1üX End V+1üV:S+XüS:T+1üT Output(5,1," "):Output(5,1,X) round(S/V,2)üM "GEMIDDELDE" Output(5,11," "):Output(5,11,M) Output(2,1,V) Goto A
KANSEN
- 20 -
5
Beroemde problemen
Het probleem van de twee kaarten Iemand heeft twee kaarten. Een ervan (kaart A) heeft twee witte zijden, de ander (kaart B) heeft een witte en een zwarte kant. Hij pakt er aselect een van en laat je een willekeurige zijde zien. Stel dat die wit is: wat is dan de kans dat ook de achterzijde wit is? Je kunt dit experiment op twee manieren bekijken.
model I
Volgens model I zijn er twee, even waarschijnlijke mogelijkheden. Ofwel het was kaart A met de achterkant wit.
W
A
W A
Ofwel het was kaart B met de achterkant zwart. Fifty-fifty dus. De kans op wit is ½.
B A
W
Z
model II Volgens model II zijn er drie, even waarschijnlijke mogelijkheden (1 is de voorkant, 2 de achterkant): Bij kaart A: (W1,W2) en (W2,W1); bij kaart B: (W1,Z2); (Z1,W2) vervalt. De kans op wit is dus 2 3 Laat dit keer het experiment beslissen!
W1 A W2
W
W1
Z
B Z2
Het schermplaatje rechts hierboven toont het resultaat na slechts tien simulaties. Dat zegt uiteraard nog niets. Maar kunnen we na bijvoorbeeld duizend simulaties misschien een conclusie trekken?
KANSEN
- 21 -
W
KAARTEN2 0üW:0üZ:0üX:0üT While T<S "KAART A HEEFT 2 WITTE KANTEN" randInt(0,1)üA "A=0 IS KAART A, A=1 IS KAART B" randInt(0,1)üV "V=0 IS WITTE VOORK, V=1 ZWART" W+(A=0)üW Z+(A=1 and V=0)üZ W+ZüT Output(1,12,T) Output(2,12,W) Output(3,12,Z) End round(100W/T,2)üP:Output(5,1,P) Output(5,7,"PROCENT") Output(6,1,"ACHTERKANT WIT") Pause :ClrHome
Het probleem van de oppas Er zijn twee kinderen om op te passen. De oppas belt aan. Een meisje opent de deur. Hoe groot is de kans dat het andere kind ook een meisje is? Is hier model I of model II op van toepassing? Of geen van beide? Probeer het antwoord via simulatie te vinden. Maar gebruik daarvoor meer dan 11 simulaties (liefst 1000). OPPAS 0üM:0üT Lbl 0 randInt(0,1)üA randInt(0,1)üB "0=MEISJE" "M TELT DE MEISJES" "T IS TOTAAL" If A+B=2:Goto 0 T+1üT If A=0:M+1üM Output(1,1," ") 100M/TüP round(P,1)üP Output(1,1,P) Output(4,1,T) Goto 0
OPGAVE 10.
Teken het hierbij behorende blokdiagram.
KANSEN
- 22 -
Het probleem van de drie deuren Je bent gast in een spelshow. Voor je zie je drie deuren, maar achter slechts één ervan ligt een (formidabele) prijs. Jij kiest een van de drie deuren. Nu komt het: de quizmaster opent één van de twee andere deuren en laat je zien dat daarachter de prijs niet ligt. (De quizmaster weet welke deur de goede is!). Nu krijg jij de kans alsnog te wisselen, dus voor die derde dichte deur te kiezen. Je kunt ook bij je keus blijven, de deur van je eerste keus. Wat doe je: wissel je of wissel je niet? Dit probleem staat bekend als het Monty Hall probleem of het Willem Ruis probleem, genoemd naar twee grote quizmasters uit de jaren tachtig. De kranten zijn volgeschreven met commentaar van lezers, waaronder veel wiskundigen en statistici. Niemand hoeft zich te schamen als hij dit ogenschijnlijk simpele vraagstuk verkeerd aanpakt. Veel wiskundigen zijn Quizmaster hem daarin voorgegaan! We bieden de lezer -net als in het probleem met de twee kaarten- twee oplossingsmodellen aan.
Kies deur: opent:
Kies deur:
1
Wissel:
1
1
3
3
½
In model I veronderstellen we dat de ongeopend 2 overgebleven deuren gelijke kans (dus ½) hebben op succes. ½ Als je bijvoorbeeld deur 2 het eerst koos en de quizmaster 3 opent deur 1, dan maakt het niet uit of je van keuze verandert (dus deur 3 kiest), want deur 2 heeft dezelfde kans op succes. Wel of niet wisselen van deur maakt dus niet uit. Aldus oplossingsmodel I. Model II. Stel dat je weer deur 2 kiest. Er is een even grote kans ( 1 3 ) dat de prijs achter één van de drie deuren ligt. Stel ook, dat je altijd van keuze wisselt. Als je prijs achter deur 2 lag, heb je pech in dat geval want je opent straks een andere deur. Kans op pech dus 1 3 . Maar als de prijs achter deur 1 lag (kans 1 3 ) loopt het goed af want de quizmaster moest deur 3 openen waarna jij van keuze wisselde en deur 1, waarachter de prijs ligt, opende. Advies is dus niet te aarzelen maar te wisselen van mening, met winstkans 2 3 .
Kies deur:
1
Prijs achter:
Wissel:
1
Raak
1/3 2
Mis
2
Mis
1/3 3
3
Raak
Het programma DEUREN3 bootst de situatie na. Je moet echter een groot aantal keer het spel naspelen voordat het verschil tussen kans 1 2 en 1 3 duidelijk begint te worden. Uiteraard is er ook een programma dat voor je speelt, een automatisering van het DEUREN3. We hebben dat programma DEURAUTO genoemd. Probeer het maar eens uit. (Waarom is het totaal gewisseld+niet gewisseld niet 100 procent?).
KANSEN
- 23 -
DEURAUTO 0üP:0üQ:0üR:0üS Lbl A randInt(1,3)üD:randInt(1,3)üW If DøW:6-D-WüT randInt(0,1)üA If D=W Then If D=1:2+AüT If D=3:1+AüT If D=2:1+2AüT End randInt(0,1)üB If T=1:2+BüE If T=3:1+BüE If T=2:1+2BüE If D=E and E=W:P+1üP If D=E and EøW:Q+1üQ If DøE and E=W:R+1üR If DøE and EøW:S+1üS Output(2,10,R):Output(2,14,S) Output(3,10,P):Output(3,14,Q) Output(4,10,P+R):Output(4,14,Q+S Output(6,1,"WINSTPERCENTAGE:") Output(7,1,"GEWISSELD:") Output(8,1,"NIET GEW.:") If R+Sù1:Output(7,12,round(100R/(R+S),1)) If P+Qù1:Output(8,12,round(100P/(P+Q),1)) Goto A
Volledige inductie De som van de rekenkundige rij 1 + 2 + 3 + … + n is gelijk aan ½n(n+1). Er is een prachtige manier om dat te bewijzen, volgens het eeuwenoude principe van de volledige inductie. Dat werkt als volgt. Voor n=1 klopt de formule, want 1 = ½.1.2 Stel even dat de formule voor een bepaalde waarde n bewezen is, dus ga uit van: 1+2+3+…+n = ½n(n+1) ……….. (*) Als we kunnen bewijzen dat de formule ook voor de opvolgende waarde n+1 geldt, dan zijn we klaar. Immers, voor n=1 klopt hij, voor de opvolger n=2 dus ook, voor de opvolger n=3 ook, enzovoorts. Hoe bewijs je dat formule (*) ook voor de opvolger n+1 geldt? Voeg links de n+1-ste term toe aan de rij en vul rechts in n+1 voor n: Is 1+2+3+…+n + (n+1) = ½ (n+1)(n+1+1) voor elke waarde van n groter gelijk 1?? Gebruik makend van de juistheid van (*): Is ½n(n+1) + (n+1) = ½(n+1)(n+2) ?? Delen door (n+1) mag, omdat n+1¹0 is: Is ½n.1 + 1 = ½.1.(n+2) ?? Ja! Want ½n + 1 = ½n + 1 voor elke waarde van n (³1) Dus, nogmaals: voor n=1 klopt formule (*) (vul maar in), daaruit volgt dat hij voor de opvolger n=2 ook klopt en voor de opvolger daarvan n=3 ook. Hij klopt dus voor alle n. Nog een voorbeeld.
KANSEN
- 24 -
2
Voor de somrij van de kwadraten geldt: 1+4+9+…+n =
1 n(n + 1)(2n + 1) 6
Bewijs: voor n=1 staat er 1 =
1 6
(1 + 1)( 2 + 1) = 1 en dat is juist. 2
Voeg links de n+1-ste term (n+1) toe en vervang rechts n door n+1: 2
2
1+4+9+…+n + (n+1) = en onderzoek of dit klopt: 1 n(n + 1)(2n + 1) + (n + 1) 2 6
1 (n + 1)(n + 6
2)(2n + 2 + 1)
= 16 (n + 1)(n + 2)(2n + 3) ???
Deel links en rechts door (n+1), dan is de vraag of het volgende klopt: 1 n(2n + 1) + (n + 1) = 16 (n + 2)(2n + 3) ??? `6 Vermenigvuldig met 6 en werk uit: 2n 2 + n + 6n + 6 = 2n 2 + 7 n + 6 dit is juist voor alle n ³ 1
Conclusie: voor n =1 is de formule juist en voor iedere volgende waarde van n ook.
Pak twee getallen uit n. Hoe groot is de grootste van de twee? (volgens de verwachting) We hebben een uitkomstenverzameling {1,2,3,…,n}. Zoals de uitkomsten van een dobbelsteen: {1,2,3,4,5,6}. Een simulatie. Pak aselect twee uit die n getallen. We kijken naar de verwachtingswaarde van het grootste getal (maximum) van de twee. Hiernaast staan in een roosterpatroon van 6 bij 6 de 7 worpen van twee dobbelstenen gemarkeerd waarvoor het maximum 4 is, namelijk (1,4) (2,4) (3,4) (4,4) (4,3) (4,2) en (4,1). Voor het trekken van 2 uit 6 met teruglegging is P(max=4) dus 7/36. 2
2
Met teruglegging is de verwachtingswaarde 1/n .S k(2k-1) = (4n +3n-1)/(6n) Deze betrekking is te bewijzen met volledige inductie. Met de GR kan het ook: sum(seq(X(2X-1),X,1,N)/N2). In het geval zonder teruglegging moet je de diagonaal (1,1) (2,2) (3,3) … weglaten. Zonder teruglegging is de verwachtingswaarde dus: 1/(n(n-1)).S 2k(k-1) = (2n+2)/3 oftewel sum(seq(2X(X-1),X,1,N)/(N(N-1))). Hierna doen we 400 simulaties met twee dobbelstenen. Dat is dus een trekking MET terugleg.
KANSEN
- 25 -
MAX2UITN Met teruglegging: Input "SIMULATIES? ",S For(X,1,S) Output(5,1,X) randInt(1,N)üA randInt(1,N)üB max(A,B)üM 1+L1(M)üL1(M) End
Zonder teruglegging: For(X,1,S) Output(5,1,X) randInt(1,N)üA Repeat AøB randInt(1,N)üB End max(A,B)üM 1+L1(M)üL1(M) End
De verwachtingswaarde (Expectation Value) voor de hoogste waarde is dus 4n 2 + 3n - 1 bij trekking van 2 uit n met terugleg 6n 2n + 2 Exp zonder ( MAX ) = bij trekking van 2 uit n zonder terugleg 3 Expmet ( MAX ) =
Welke formules krijg je voor de laagste waarde van de twee? Merk eerst het volgende op over de verschillende kansen: P(MAX=k) = P(MIN=n+1-k), hier geïllustreerd voor n=6 en k=3. hetgeen in beide gevallen leidt tot de conclusie: Exp( MIN ) + Exp( MAX ) = n + 1
waarna het antwoord op de vraag gegeven kan worden. 2n 2 + 3n + 1 6n n +1 Exp zonder ( MIN ) = 3 Expmet ( MIN ) =
KANSEN
- 26 -
OPGAVEN 11.
a) Uit het roosterpatroon hiernaast kun je de formule 2
1+3+5+7+…+(2n-1) = n aflezen. Hoe? Bewijs deze formule ook met volledige inductie. Aanwijzing: de volgende term is 2n+1.
b) En hoe zie je aan het roosterpatroon in het tweede plaatje dat de somformule voor de rekenkundige rij 1+2+3+...+n s n = 1/2n(n+1) deugt?
c) Verander in het programma max(A,B) in min(A,B) en controleer experimenteel de formules voor de verwachtingswaarde van de laagste trekking. 12.
Herhaal het getoonde experiment voor trekking van 2 uit 6 zonder terugleg.
13.
Bereken met kansrekening de kansen (met/zonder) voor een trekking van 2 uit 4 en doe een paar honderd simulaties ter controle.
14.
Bewijs de betrekking 1/(n(n-1)).S 2k(k-1) = (2n+2)/3 met volledige inductie.
Eenentwintigen De kans om een aantal ogen precies te gooien. Bij het dobbelspelletje Makao probeer je zo dicht mogelijk bij totaal 9 ogen te komen. Wie de negen overschrijdt valt af. Wie precies op negen uitkomt heeft altijd gewonnen. Hiermee verwant zijn de spellen Eenentwintigen en Black Jack. We simuleren in TREFKANS het volgende. Je gooit telkens met een dobbelsteen, telt op, en probeert met de totaalscore precies op een doelwaarde N uit te komen. Hoe groot is de kans P(N=n) dat dit lukt? (Bij overschrijding van N moet je overnieuw beginnen). Voor de doelwaarde 1 is die kans P(N=1)=1/6. Voor N=2 zijn er twee wegen die naar succes leiden, namelijk: direct 2 gooien met kans p(2)=1/6 of indirect via twee keer een 1 met p(1,1)=1/36 De totale trefkans P(N=2) is dus 1/6+1/36=7/36=0,19444 En P(N=3) is p(3)+p(1,2)+p(2,1)+p(1,1,1)=49/216=0,22685 Bij P(N=4) hebben we een paar grafen getekend die je van links naar rechts moet lezen. Plaatje 1. Je kunt precies 4 gooien in één keer. De kans daarop is p(4)=1/6. Plaatje 2. Je kunt via 1 en daarna 3 met kans p(1,3) = 1/36 ook op 4 uitkomen. Totaal 1 plus 2: 1/6 + 1/6´P(N=1) Plaatje 3. Via N=2 is het p(2,2) + p(1,1,2) Totaal: 1/6 + 1/6´P(N=1) + 1/6´p(N=2) Plaatje 4. Via N=3: p(3,1) +p (1,2,1) + p(2,1,1) + p(1,1,1,1) Totaal: zie opgave 32.
KANSEN
- 27 -
plaatje 1
1
plaatje 2
2
3
4
1
3
4
2
3
4
plaatje 4
plaatje 3
1
2
2
3
4
1
Voor grote doelwaarden is het niet gemakkelijk een formule te geven. Wel is duidelijk dat P(N=n) voor grote waarden van n naar een limiet convergeert (ongeveer 0.286). In opgave 27 wordt gevraagd deze limiet exact (dat wil zeggen als niet-decimale breuk) te berekenen. Het programma TREFKANS simuleert 1000 keer en berekent (recursief) de theoretische trefkansen t/m P(N=50). OPGAVEN 15.
Geef een recursieve formule die voor n£6 P(N=n) uitdrukt in P(N=1) t/m P(N=n-1) en voor n³7 uitdrukt in P(N=n-6) t/m P(N=n-1).
16.
Geef voor n£6 een directe formule die P(N=n) uitdrukt in n. Zie ook de grafiek hiernaast.
17.
Bepaal de exacte limietwaarde voor grote N. Gebruik daarvoor het eerder (bij het programma COMPLEET) behandelde principe: verwachting =
TREFKANS For(T,1,1000) While S
n
åc × k = c × åk k =1
k =1
Daar gaan we. d = [0 × 1 n + 1 × 2 × n -1 n + 2 × 2 × n - 2 n + 3 × 2 × n - 3 n + ...] = n2 [1 × ( n - 1) + 2 × ( n - 2) + 3 × (n - 3) + ...] n -1
n -1
(3)
n -1
= n2 ´ å k × (n - k ) = 2n å k - 2 å k 2 k =1
k =1
(1) en( 2)
=
2 × 1 2 (n - 1)n - 2 × 1 6 (n - 1)(2n - 2 + 1)
k =1
Tenslotte komt er: d = n(n - 1) - 1 3 (n - 1)(2n - 1) = (n - 1)[n - 1 3 (2n - 1)] = (n - 1)( 1 3 n + 1 3 ) = 1
3 ( n - 1)( n + 1)
=
1
3 (n
2
- 1)
Zodat de formule voor de verwachte gemiddelde afstand van een rij voor en na het grondig schudden gevonden is: d =
1
3 (n
2
- 1)
We controleren het behandelde geval met n =4 nog even: d = kunnen zijn.
1
3 (16 - 1)
= 5 waarmee we gelukkig
Dan wordt het nu tijd om de achterliggende theorie te toetsen met de praktijk. Simuleren dus maar. De perfectionist die op zoek gaat naar de maximumafstand tussen L1 en L2 zal ontdekken dat perfecte wanorde soms gelijk is aan perfecte orde! Dus zover zullen wij maar niet gaan. Voor n = even namelijk is d MAX = 1 2 n 2 ; voor n = oneven is d MAX = 1 2 (n 2 - 1) . Die "ideale" schudsituatie hoort bijvoorbeeld bij de omgekeerde rij {n, ... ,3,2,1}; het schudden is dan 2 te lang doorgegaan. Ik laat het aan de lezer over, om zelf een norm te bepalen voor d tussen 0×(n -1) 2 en ½×(n -1) die de mate van het door elkaar schudden vastlegt. In het programma SCHUDRAN wordt een randomrij met lengte N geproduceerd waarvan de nog niet gebruikte getallen in L2 worden opgeslagen, net zo als dat eerder bij het trekken zonder terugleg gebeurde. De matches komen in L3 en de afstand in A.
KANSEN
- 32 -
Bij de uitvoering van SCHUDRAN zal blijken dat deze manier van "schudden zonder terugleg" bij lange rijen (n groter dan ongeveer 15) enorm traag verloopt. Het voordeel van de randomschudding is echter dat het gesimuleerde aantal matches en de gemiddelde afstand aardig voldoen aan de theoretische voorspellingen. Bij N=10 en slechts 182 simulaties kloppen de experimentele waarden voor m en d al erg aardig met de theoretische waarden 1 en 33. Zie het schermplaatje. Lbl HH Fill(0,L2) 0üT:H+1üH Lbl 0 randInt(1,N)üR If sum(L2=R)=0 Then T+1üT:RüL2(T) End If T
Het tweede schudprogramma SCHUDIN verwisselt herhaaldelijk (aselect) twee getallen. Ongeveer zoals een onhandige kaartspeler dat zou doen.
Lbl 0 W+1üW randInt(1,N)üX randInt(1,N)üY L2(X)üZ L2(Y)üL2(X) ZüL2(Y) sum(abs(L1-L2))üA sum(L1=L2)üM Output(1,10,M) Output(2,10,A) Output(8,10,W) If A
"KIES 2 WILLEKEURIGE POSITIES" "EN VERWISSEL DE DAARBIJ" "BEHORENDE GETALLEN IN L2" "A=AFSTAND"
Het schudden "kaart-voor-kaart" via de verwisseling van telkens twee "kaarten" gaat met een constante snelheid. Aanvankelijk (bij kleine waarden van N) vinden we dat nogal traag. Maar bij grote waarden van N (bij een spel van 52 kaarten bijvoorbeeld) verdient zich dat terug. In de laatste programmaregel hebben we de norm voor voor een aanvaarbare mate van schudding gelegd bij een afstand van 1 3 n 2 ; je zou ook met minder wanorde genoegen kunnen nemen (bijvoorbeeld bij een 2 afstand ¼ n ) of je zou kunnen kijken naar het aantal matches.
KANSEN
- 33 -
Matches Simuleer het schudden van N getallen in een niet al te lange rij {1,2,...,N} en kijk naar het aantal matches (niet de lucifers). De theorie is enkele bladzijden hiervoor beschreven. We hebben het gemiddeld aantal matches meteen maar even uitgerekend (onder THEORIE).
MATCHES 0üT For(U,1,S) Output(1,5,U) Lbl 0 randInt(1,N)üR If sum(L2=R)=0 Then T+1üT:RüL2(T) End If T
KANSEN
"SCHUDDEN" "VOLGENS" "SCHUDRAN"
- 34 -
6
Continue kansvariabelen
De kanssimulaties tot dusver betroffen voornamelijk discrete situaties, waarbij de waarnemingsgetallen (X) gehele waarden aannamen, zoals het aantal ogen bij een dobbelsteen. Ingewikkelder wordt het bij continue kansexperimenten, met x-waarden die niet per definitie geheel zijn. Ik gebruik hier met opzet de cursieve, kleine letter x als continue variabele.
Een duel (DUELPQ) A en B houden een duel en schieten om de beurt A schiet met een trefkans P; B met een trefkans Q. A begint. Wie het eerst raak schiet, heeft gewonnen. Wat is de kans dat A wint? ClrHome Prompt P,Q 0üT:0üU Lbl 0 randüA If A
De theorie achter deze simulatie berust op een meetkundige rij. De kans op winst van A wordt immers berekend via het patroon RAAK of MIS-MIS-RAAK of MIS-MIS-MIS-MIS-RAAK of ... Dus: 2 2 3 3 p + (1-p)(1-q)p + (1-p) (1-q) p + (1-p) (1-q) p+ ... = (volgens de somformule van een oneindige meetkundige rij) 1 p p q = . = 1 + 1 -1 1 - (1 - p )(1 - q ) p + q - pq q p Duidelijk is dat A hierbij het voordeel van de eerste zet heeft. De winstkansen voor A en B zijn gelijk als de laatste breuk gelijk aan ½ gesteld wordt. Dit leidt tot: 1 - 1 = 1. p q Probeer maar eens met bijvoorbeeld p = 0.20 en q = 0.25.
Een meting van p, met darts We kunnen de waarde van p (pi) op een natuurkundige manier "meten". De oppervlakte van een figuur is evenredig met het kwadraat van een lengte (of straal, in dit geval). Ga dus uit van de oppervlakteformule voor een cirkel: 2 p.r . We werpen een groot aantal keren met een denkbeeldig pijltje (dart) naar een 2 eenheidscirkel, dat is een cirkel om O(0,0) met straal 1 (en oppervlakte p.1 = p). Om die eenheidscirkel is een vierkant beschreven (het dartsbord) met zijde 2 en oppervlakte 4, zodat de kans (bij lukraak gooien) om binnen de cirkel uit te komen p/4 is.
KANSEN
- 35 -
Met -1+2rand krijg je randomgetallen tussen -1 en 1. We noteren een treffer, als de afstand tot de 2 2 oorsprong (0,0) kleiner dan 1 is, dus x +y <1. Je zult merken dat deze meting nogal onnauwkeurig is, zelfs als je 5000 keer simuleert. SIMULPI 0üT "T=TELLER" For(I,1,N) ú1+2randüX "RANDOM TUSSEN -1 EN 1" ú1+2randüY If XÜ+YÜ÷1 "PIJLTJE BINNEN CIRKEL" Then T+1üT "TELLEN" End Output(6,11,I):Output(7,11,T) round(4T/I,3)üP Output(8,4,P) End
OPGAVE 19. Schrijf een programma dat de darts treffers op de hiernaast getekende wijze illustreert na 2000 worpen. De stippen moet je achter elkaar zien verschijnen. Het duurt wel een paar minuten.
Ontwerp zelf een randomgenerator Een randomgenerator is een eindeloos voortgaande rij willekeurige getallen tussen 0 en 1, recursief opgebouwd volgens het principe:
ì x0 = kiem í î xn +1 = f ( xn ) Voor het startgetal kiem nemen we een getal dat niet te dicht bij nul ligt, bijvoorbeeld Ö11. Voor de hand liggend is om voor f een breukdeel (Fraction-Part) te gebruiken omdat het decimale gedeelte van een getal altijd tussen 0 en 1 ligt. Nou zijn deze getallen xn bepaald niet aselect ontstaan. Elk getal volgt volgens een vast voorschrift uit zijn voorganger. Afhankelijker dan dat kan het haast niet. Vandaar dat we het woord pseudo-random gebruiken. Ik doe maar eens wat, f (x)= fPart(1234×p× x 1/7) bijvoorbeeld. Deze doe-het-zelf randomgenerator noem ik een pseudo-pseudo-randomgenerator, omdat hij gemaakt is door een beunhaas (schrijver dezes). Van de uitkomsten kun je een klassenverdeling maken (de x is hieronder in tien klassen verdeeld, 0-0.1 t/m 0.9-1.0) om na te gaan of alle getallen tussen 0 en 1 wel ongeveer even vaak aan de beurt komen. En dat is aardig gelukt, kijk maar. Input "KIEM?",X For(N,1,2000) fPart(1234ÄX^(1/7))üX iPart(10X)üI L2(I+1)+1üL2(I+1) Output(8,1,N) End
KANSEN
- 36 -
De professionele pseudo-randomgenerator rand zorgt voor willekeurige (aselecte) getallen tussen 0 en 1. Kies eerst een startwaarde (nul bijvoorbeeld) met 0®rand en druk een paar keer op ENTER. Zo krijg je op alle modellen TI-83 en TI-84 dezelfde reeks randomgetallen (zie het schermplaatje), omdat al deze modellen dezelfde microprocesser (Z-80) ingebouwd hebben. Bij een andere startwaarde wordt een ander rijtje "gegenereerd". Andere microprocessors, andere computers en andere startwaarden leveren andere randomgetallen. Dat maakt niet uit, mits de geproduceerde getallen -hoe zal ik het zeggen- een beetje normaal verdeeld zijn. Normaal verdeeld is eigenlijk een buitengewoon onhandige woordkeuze, de suggestie wekkend dat de getallen in het midden (rond 0.5) vaker voorkomen dan aan de rand. Een beetje regelmatig verdeeld, is dat een betere formulering? Nee, want hierbij denk je wellicht dat de onderlinge verschillen min of meer constant zijn, het rijtje 0.28 0.38 0.48 0.58 (ik neem voor het gemak maar twee decimalen) wekt bepaald niet een indruk van willekeurigheid. Afwisselend misschien? Maar de afwisselende getallen 0.28 0.72 0.23 0.77 0.21 0.79 zijn duidelijk niet aselect ontstaan! Tel ze maar eens op. Het is moeilijk een geschikte term te vinden voor de bedoelde aselectheid. Wiskundigen, liefhebbers van ondubbelzinnigheid, gebruiken de woorden homogeen of uniform. Een goede lukraakrij moet op zijn minst voldoen aan een vooraf te berekenen gemiddelde en een bepaalde spreiding (variantie). De zogenaamde kansdichtheid van de randomfunctie rand is gedefinieerd door f (x) = 1 voor x tussen 0 en 1 en f(x) = 0 voor andere x-waarden.
Het woord kansdichtheid is ook weer niet exact, het suggereert immers dat elk randomgetal een kans (ongelijk nul) heeft. Maar de kansen op 0.9435974025 en 0.908318861 zijn beide 1 op de10^10 (de rekennauwkeurigheid van het machientje is 10 decimalen) en je kunt hier evenmin van een kans spreken als bij "de kans dat iemand precies 1 meter 83,76 cm lang ben is ...". Een lijn heeft immers geen breedte. Wel zinvol is te kijken naar x-waarden tussen twee grenzen, een intervalletje Dx. De oppervlakte van een staafje met deze breedte is f(x).Dx = 1.Dx; de totale oppervlakte van alle staafjes 1
in het eenheidsvierkant is S(1.Dx)=1, in limiet hebben we daar de integraalnotatie ò 1. dx = 1 voor. 0
De kans P(X
KANSEN
- 37 -
Met het begrip kansverdeling (distribution) valt beter te werken. De kansverdelingsfunctie, hier aangeduid met de hoofdletter P, R
ò f ( x) dx ,
P( X < R) =
-¥
is te vergelijken met de cumulatieve discrete kansverdeling k
P ( X £ k ) = å pi . i =1
De verwachtingswaarde E van X (E komt van expectation. m is de Griekse letter m), +¥
m = E( X ) =
ò x. f ( x) dx -¥
is te vergelijken met het gemiddelde in het discrete geval X =
n 1 n xi . f i of X = å xi . pi . å n i =1 i =1
2
De verwachtingswaarde van X is dus +¥
E( X 2 ) =
òx
2
. f ( x) dx .
-¥
2
En de variantie var(X) = s (X) als spreidingsmaat, volgend uit var ( X ) = E ( X - E ( X )) 2 , is analoog aan het discrete geval met 2
2
2
s 2 = ( X - X )2 = ( X 2 - 2 X X + X ) = X 2 - 2 X ×X + X = X 2 - X . 2
De standaarddeviatie s is uiteraard de wortel hieruit: s = X 2 - X . PANORAMA R
P( X < R) =
ò f ( x) dx ; de kansverdeling is een integraal van de kansdichtheid -¥
dus de kansdichtheid is de afgeleide van de kansverdelingsfunctie 1
+¥
x=
]
1
2 ò x. f ( x) dx dus R = x.1dx = 1 2 x 0 = 1 2 ; m = 0.5
ò
-¥
0
1
R 2 = ò x 2 .1dx = 13 x3
]
1 0
=
1
3
en
0
2
2
s = R2 - R =
KANSEN
1
3
- 14 =
1 ; 12
s 0.2887
- 38 -
Ik ben nieuwsgierig naar de kwaliteit van mijn pseudo-pseudo randomgenerator, vergeleken met de professionele randommaker. En dat valt, aan het histogram te zien, best mee! Het professionele plaatje (pseudo generator) lijkt zelfs wat minder mooi omdat rondom het midden een beetje een deuk zit. Het oog wil ook wat, natuurlijk, maar op getallen en berekeningen kun je blindelings vertrouwen, om een bij het onderzoek passend woord te gebruiken. We gaan dus het gemiddelde en de variantie berekenen na 4 keer 500 simulaties van de twee beschouwde randomgeneratoren en die vergelijken met de theoretische waarden. theorie m = 0.5 s = .289
pseudo random m = .490 s = .293 m = .497 s = .296 m = .498 s = .294 m = .515 s = .281 na 2000 simulaties: m = .500 s = .291
pseudo-pseudo random m = .480 s = .295 m = .509 s = .286 m = .497 s = .288 m = .500 s = .284 na 2000 simulaties: m = .497 s = .288
Nou blijven we nog met één hoogst fundamentele vraag zitten. De randomgetallen die volgens het recursieve model zijn geproduceerd, zijn per definitie afhankelijk. En onafhankelijkheid is beslist een noodzakelijke voorwaarde voor verantwoorde kansexperimenten. Onafhankelijkheid kan getoetst worden met de productregel. Als voor opvolgende getallen op den duur (na een groot aantal simulaties) geldt: rand ´ rand = rand ´ rand dus rand ´ rand = 12 ´ 12 = 14 dan zullen we aannemen dat de randomgenerator ook in dit aspect bevredigend werkt. En kijk:
leverde na een tijdje:
Ook na 5000 tests kwam ook de pseudo-pseudo-randomgenerator tot de hoogst bevredigende waarde 0.2514 voor het gemiddeld product van opvolgende randomgetallen. Waarmee dit onderzoek naar de betrouwbaarheid van randomgenerators afgesloten wordt. Laten we gauw beginnen te experimenteren met continue kansvariabelen.
De afstand tot een cirkelrand Punt P ligt willekeurig binnen een cirkel met straal 1. Wat is de verwachtingswaarde van de (gemiddelde) afstand van P tot de cirkelomtrek? Programmanaam: PINCIRKL.
KANSEN
- 39 -
0üT:0üN Lbl 0 2rand-1üX:2rand-1üY ð(XÜ+YÜ)üD If D>1:Goto 0 1-DüA:N+1üN:T+AüT Disp T/N Goto 0
De theoretische achtergrond. Stel a is de afstand van P tot de cirkel. De kansverdeling wordt gegeven door een oppervlakte verhouding: opp ring p ×12 - p × (1 - a ) 2 = = 2a - a 2 2 opp eenh cirkel p ×1 De kansdichtheid is de afgeleide hiervan, namelijk de gemiddelde afstand tot de rand van de cirkel is 1
a = ò a × (2 - 2a ) da = a 2 - 23 a 3 0
]
1 0
d da
a 1-a
(2a - a 2 ) = 2 - 2a dus
= 13
OPGAVE 20.
Bereken de kans dat een randompunt binnen de eenheidscirkel op een afstand kleiner dan ½ ligt van de omtrek. Controleer met simulatie of deze berekening klopt.
De afstand tot de rand van een vierkant (AFP4KANT) P ligt random binnen een vierkant met zijde 1. Wat is verwachtingswaarde van de (kortste) afstand van P tot de zijden van dat vierkant? In het programma veronderstellen we dat het vierkant in het eerste kwadrant ligt met O als een van de hoekpunten en dat P de coördinaten (x,y) heeft. ClrHome 0üS:0üT Lbl 0 randüX:randüY min({X,Y,1-X,1-Y})üD T+1üT: S+DüS Disp S/T Goto 0 1 a P
1-y 1
a x
1-x
P
1 -2a 1
y 0
KANSEN
1
- 40 -
Stel a is de afstand van P tot een zijde. Merk op dat a £ ½ is. P ligt tussen de twee vierkanten in. Het 2 2 binnenste vierkant heeft een zijde 1-2a dus een oppervlakte (1-2a) = 1-4a+4a . De kans dat P op een 1 - (1 - 4a + 4a 2 ) afstand kleiner dan a van de omtrek ligt is = 4a - 4a 2 . De kansdichtheid is de 1 1
2
afgeleide hiervan, dus 4-8a. De gemiddelde afstand is dus a = ò a × (4 - 8a ) da =2a 2 - 83 a 3 0
]
1
2
0
= 16 .
OPGAVEN 21.
Bereken de kans dat een randompunt binnen het eenheidsvierkant op een afstand kleiner dan 0.2 van de rand ligt. Controleer met simulatie of deze berekening klopt.
22.
Een punt beweegt random binnen een rechthoek van 3 bij 2. Noem a de afstand van dat punt tot de dichtstbijzijnde zijde van die rechthoek. a) bepaal (door theoretische berekening en met simulatie) de kans dat a < 0.5 is. b) laat zien (theoretisch en via simulatie) dat de gemiddelde afstand van P tot de dichtstbijzijnde zijde 7/18 is.
a a 2-2a
2
3-2a
P
3
De stomphoekige driehoek Binnen de eenheidscirkel (middelpunt O, straal 1) liggen, willekeurig, twee punten A en B. De vraag is, hoe groot de kans is dat driehoek AOB stomphoekig is. Geef de punten de coördinaten A(P,Q), B(R,S) en O(0,0). De snelste manier om te constateren of een 2 2 2 hoek stomp is is met het inproduct. Niet iedereen weet dat, maar met de regel c >a +b kan het ook. 2
2
2
2
2
2
2
2
(1) Als O een stompe hoek is geldt AB > OA + OB dus (R-P) +(S-Q) > P +Q +R +S 2 2 2 2 2 2 2 2 2 (2) Als B stomp is geldt OA > OB + AB dus P +Q > R +S +(R-P) +(S-Q) 2 2 2 2 2 2 2 2 2 (3) Als A stomp is geldt OB > OA + AB dus R +S > P +Q +(R-P) +(S-Q)
2
ClrHome 0üT:0üU Lbl 0 U+1üU 2rand-1üP:2rand-1üQ:2rand-1üR:2rand-1üS PÜ+QÜüB:RÜ+SÜüA (R-P)Ü+(S-Q)ÜüC If B>1 or A>1:Goto 0 If C>A+B or B>A+C or A>B+C: T+1üT Disp T/U Goto 0
Het vermoeden rijst dat het antwoord ¾ is. Econometriestudent Tran Phi Hung van de Vrije Universiteit in Amsterdam kwam met het volgende wonderschone (dus korte) bewijs dat dit vermoeden juist is. KANSEN
- 41 -
Noem het verst van O liggende hoekpunt A en bekijk de cirkel C1 met middelpunt O die door A gaat. Hoek A kan nu niet stomp zijn. Er zijn twee gebieden waar het andere hoekpunt (B) kan liggen, ingeval de driehoek stomphoekig is en OB
Dating Twee geliefden proberen elkaar dagelijks te ontmoeten voor een innige omhelzing. Beiden arriveren random tussen 17 en 18 uur bij het meeting point. Afgesproken is om hoogstens 10 minuten te wachten op de ander. Hoe groot is de kans op een omhelzing?
Stel de aankomsttijden zijn resp x en y minuten na 17.00 uur. Bekijk nu het kansvierkant met 2 oppervlakte 60´60 = 3600 min (vierkante minuten dus). Elk punt (x, y) daarbinnen staat voor één van de mogelijke gebeurtenissen. Het tijdsverschil tussen de geliefden x-y of y-x moet (liefst) kleiner zijn dan 10 minuten: voor het ontmoetingsgebied geldt dus |x-y|<10. De lijn y = x beschrijft gelijktijdigheid; het gebied buiten de lijnen y = x + 10 en y = x – 10, de twee rechthoekige driehoeken, 2 bestrijkt een oppervlakte van in totaal 2500 min . 2 Voor de gelukzalige momenten zijn dus 3600-2500=1100 min over. Met andere woorden: de kans op 1100 11 een omhelzing is = » 0,3056 . Deze laatste eis wordt in tibasic weergegeven met vooraf 0üS 3600 36 en daarna: S+(abs(60rand-60rand)<10)üS Een vergelijkbaar vraagje is: Trek twee keer een randomgetal tussen 0 en 1. Hoe groot is de kans dat het product minder is dan 0,5?
Aangezien x× y < ½ is, moeten we binnen het kansvierkant de oppervlakte onder de kromme 0,5/x in rekening brengen:
KANSEN
1 1 1 + dx = 2 0ò,5 2 x
1
2 (1 + ln 2)
» 0,8466 en dat is de gevraagde kans.
- 42 -
DATING 0üS:0üT Lbl 0 T+1üT S+(abs(60rand-60rand)<10)üS Output(1,11,T) Output(4,10,round(S/T,3) Goto 0
RANDPROD 0üS:0üT Lbl 0 T+1üT S+(rand*rand<.5)üS Output(1,11,T) Output(4,10,round(S/T,3) Goto 0
Maak een driehoek van een stok Martin Gardner beschrijft in "The Colossal Book of Mathematics" in hoofdstuk 21 en 22 hoe verraderlijk sommige ogenschijnlijk simpele kansberekeningen zijn. Befaamde wiskundigen als Leibniz, d'Alembert en Erdös zijn ons voorgegaan in het maken van fouten. Eerstgenoemde dacht dat de kans op 12 even groot is als de kans op 11 bij het werpen met twee dobbelstenen; d'Alembert veronderstelde dat na het werpen van erg veel keer "Munt" achter elkaar de kans op "Kop" groter zou zijn, leidend tot de conclusie dat een dobbelsteen een geheugen heeft. Een bekend probleem is dat van de stok die at random (aselect) in drie stukken geknipt wordt. Hoe groot is de kans dat uit die drie stukken een driehoek gevormd kan worden? [De driehoeksongelijkheid zegt, dat de langste zijde van een driehoek altijd korter is dan de som van de andere twee zijden]. Stel dat onze stok 1000 millimeter is. Hij wordt doorgeknipt op twee plaatsen. Methode 1: knip tegelijkertijd op twee willekeurige plekken .Methode 2: knip de eerste keer willekeurig, pak blindelings (willekeurig) een van de twee stukken, en knip een tweede keer. Methode 3: hetzelfde als methode 2, echter pak voor de tweede knip het grootste stuk van de twee; dat geeft een grotere kans maar gaat natuurlijk ten koste van de aselectheid (je selecteert het grootste stuk, omdat je daarmee twee keer zoveel kans krijgt op succes). Het programma STOKIN3 biedt methode 1 en 3 aan (de volledig aselecte experimenten) in een stapvoor-stap demo-versie en een doorlopende versie, die we echt wel nodig hebben want pas na vele honderden simulaties komt er enige lijn in. Alle beeldschermwisselingen kosten erg veel rekentijd, vandaar dat het "kale" programma STOKSNEL gemaakt is voor degene die graag zoveel mogelijk simulaties wil doen (en over een SilverEditie beschikt waarmee ruim twee keer sneller gerekend kan worden).
KANSEN
- 43 -
De theorie bij STOKIN3 (1) Twee keer tegelijk knippen. Knip tegelijk (random, uiteraard) op twee plaatsen: x en y (neem even y < x). Dat geeft een punt (x, y) in het kansvierkant met de volgende eigenschappen: x > ½ en x – y < ½ en y < ½ Het hierbij behorende gebied is de gearceerde driehoek met een oppervlakte die eenachtste is van het eenheidsvierkant. Wanneer je de tweede mogelijkheid (y > x) erbij neemt geeft dat nog zo'n driehoek zodat de kans dat van een stok, die tegelijk op twee plaatsen gebroken wordt, een driehoek gemaakt kan worden gelijk is aan ¼. (2) Knip een keer, daarna in het grootste stuk nog een keer. Knip eerst op een willekeurige plek x (stel eerst weer even dat x groter is dan ½). Knip het grootste stuk (tussen 0 en x) willekeurig in twee stukken. Er geldt: x – ½ < y < ½ . De lengte van het toegestane gebied voor y is ½ – (x – ½) = 1 – x en omdat y random ligt tussen 0 en x geldt: y = x*rand dus y × x < 1 – x oftewel y < 1-x x De oppervlakte van het toegestane gebied is dan 1
ò 1
2
1 1- x x
dx = ò ( 1x - 1) dx = - 1 2 + ln 2 » 0,19315 1
2
In het geval x < ½, is het toegestane gebied x + ½ - ½ = x. En dan zal y = (1 – x)*rand moeten zijn dus y ×(1-x) < x oftewel y < 1-x x met 1
2
ò 1-xx dx met hetzelfde
de voorwaarde ½ < y < x + ½ en als oppervlakte
0
antwoord. In totaal is de kans op een driehoek dus 2 ´ ( - 1 2 + ln 2) = -1 + 2 ln 2 » 0,3863 .
KANSEN
- 44 -
De knekelman Een stok wordt in tweeën gebroken. Wat is de verwachtingswaarde van het quotiënt
kleinste stuk ? grootste stuk
STOKIN2 ClrHome 0üT:0üQ Lbl 0 T+1üT randüR If R>.5:1-RüR Q+R/(1-R)üQ Disp Q/T Goto 0 Als de stok een lengte 1 heeft (waarom niet?) en k de lengte van het kortste stuk is, gaat het hier om het quotiënt k / (1-k). 1
Voor 0 < x < ½ moet dus de oppervlakte
2
x
ò 1 - x dx berekend worden; voor 0
½ £ x £ 1 krijg je een andere integraal, maar dezelfde oppervlakte: 1 1 (1 - x) 1 dx = ò x ò ( x - 1)dx . 1 1 2
2
De totale oppervlakte, en daarmee de verwachtingswaarde van het quotiënt
kleinste , wordt dus grootste
2 × [ln x - x ]11 = -1 + 2 ln 2 » 0,3863 . 2
Een aardige, simpele oplossing dus. Maar wiskundigen weten bij iedere aardige, simpele oplossing wel weer een nieuw probleem te bedenken. Ik citeer uit The Mathematics of Oz, Mental Gynastics from Beyond the Edge van Clifford Pickover (Cambridge U.P., 2002): The Bone Man Problem Dorothy and Dr. Oz peer into a deep hole in the ground. The bone man comes closer and opens and closes his mouth spasmodically. "In the pit", he says, "are 10 000 leg bones. I have cracked each bone at random into two pieces by throwing them against a rock. What do you think is the average ratio of the length of the long piece to the length of the short piece for each time I crack a bone? You can reason from a purely theoretical standpoint. If you cannot find the solution within two days, I will add your leg bone to the pit." Let op de subtiele omkering van het quotiënt. In het geval kleinste/grootste hebben extreem kleine stukjes vrijwel geen invloed op het verwachte quotiënt (die staan namelijk in de teller). Maar in het probleem van de Bone Man is de breuk omgekeerd en hebben extreem kleine stukjes (in de noemer) een extreem grote invloed op de waarde van de breuk! En dan hebben we echt een probleem.
KANSEN
- 45 -
Bekijk het bovenste plaatje na 100 000 simulaties voor de verhouding langste/kortste van de Bone Man. Horizontaal het aantal simulaties (n), verticaal het verloop van het voortschrijdend gemiddeld quotiënt Q . Het lijkt erop dat Q erg graag tot rust wil komen tussen de 20 en 30 maar op onverwachte momenten telkens weer opgeblazen wordt (door een noemer die akelig dicht bij nul zit). In mijn experiment moest Q zich vijf keer herstellen na een explosie. Rond simulatie nummer 80000 was er zelfs een (of meer, dat is niet te zien) soort cambrische explosie, tot buiten het schermbereik, vermoedelijk ver boven de 100. Aan de puntenwolk van het volgende plaatje is te zien dat Q nogal schommelt (n is horizontaal op een logaritmische schaal gezet). Zo registreerde ik na twee keer 10 miljoen simulaties de eerste keer Q = 53.3 en de tweede keer Q = 32.2 (de twee plusjes uiterst rechts) en na een bepaalde serie van 10 000 de uitschieter Q = 87.8 . Logaritmische regressie via [STAT] 9:LnReg leverde als best passende formule op Q = 4.4 + 2 ln n ; zie de rechte lijn. De breuk kortste/langste streeft naar een limiet van 0,3863 zoals we zagen. Waarom streeft het omgekeerde van die breuk langste/kortste niet gewoon naar het omgekeerde daarvan: 1 » 2.59 ??? 0.3863 Om daar iets van te begrijpen schakelen we terug naar het begrip verwachtingswaarde. Daarvoor moeten de Q-waarden vermenigvuldigd worden met hun kans. De randomwaarden X liggen tussen 0 en 1. Stel dat we een verdeling hebben in tien klassen, met de klassenmiddens 0.05 t/m 0.95. Alle randomwaarden (X) in de eerste klasse, tussen 0 en 0.1 vervangen we door één enkele waarde, 0.95 m=0.05; het quotiënt hierbij is Q = = 19 . De waarschijnlijkheid van zo'n quotiënt is 0.1 (er zijn 0.05 immers tien even waarschijnlijke mogelijkheden). Na n simulaties verwachten we in het ideale geval n klassen – in elke klasse één random uitkomst. Zo ontstaat voor n=10 de volgende tabel voor de berekening van de verwachtingswaarde. X 0-.1 .1-.2 .2-.3 .3-.4 .4-.5
m 0.05 0.15 0.25 0.35 0.45
Q 19 5.67 3 1.86 1.22
p 0.1 0.1 0.1 0.1 0.1
Q.p 1.9 0.567 0.3 0.186 0.122 ------- + 3.075
dus Q = 2 ´ 3.075 = 6.15 oftewel Q = (2/N)*sum(seq((1-X)/X,X,.5/N,.5,1/N)). Voor n=100 komt er Q =10.75 en voor n=999 Q =15.36.
KANSEN
- 46 -
De som S Q.p moet, als we de indeling in n klassen verder verfijnen, omgezet worden in de integraal: 1
2×
2
ò
1
2n
[
1
2 1 1- x 1 dx = 2 × ò ( - 1) dx = 2 × [ln x - x ]12 x x 2n 1
= 2 × ln
2n
1
1
2 - 2 - ln
1 2n
]
[
+ 21n = 2 × - ln 2 - 1 2 + ln 2n + 21n
]
1 1 = -1 - 2 ln 2 + 2 ln 2 + 2 ln n + n n » -1 + 2 ln n , als het aantal klassen/simulaties n ® ¥ gaat. = -2 ln 2 - 1 + 2 ln 2n +
Het is dus onmogelijk om de verwachtingswaarde te berekenen! Het venijn zit 'm in de gearceerde staart, lim ln x gaat namelijk op een x ®0
exponentiële wijze naar min-oneindig. Merk op dat de trendlijn Q = 4.4+2ln n door de puntenwolk flink wat hoger ligt dan (maar wel dezelfde helling heeft als) de verwachte waarde Q = -1+2ln n. Als je de vier hoogste Q-waarden weglaat klopt het iets beter: Q = 1.5+2ln n. Wellicht is de mediaan of het gemiddelde van de mediaan een betere maat voor het centrum van de uitkomsten. Zoals bekend, is de mediaan minder gevoelig voor extreme uitkomsten. STOK2MED ClrHome:ClrAllLists Output(1,1,"SERIENR:") Output(2,1,"MED=") For(M,1,10) of: For(M,1,100) Output(1,10,M) For(N,1,100) of: For(N,1,999) randüR If R<.5:1-RüR R/(1-R)üQ QüL1(N) End 1-Var Stats L1 MedüX:XüL2(M) Output(2,5,round(X,3)) End 1-Var Stats L2 Output(4,1,"MEDIAAN GEM :") Output(5,1,Ë)
Aldus vond ik 3.06 als gemiddelde mediaan, zowel bij n = 100 als bij n = 999. Maar daarmee is de afhankelijkheid van n verdwenen. Zo is er gelukkig weer iets om over na te denken.
Onnodige ongerustheid Een stelling die aan de stelling van Bayes voorafgaat luidt: P(B) = P(B|A).P(A) + P(B|A').P(A') De juistheid hiervan bewijzen we met enkele van de eerder gegeven stellingen. Ga maar na:
KANSEN
- 47 -
B = BÇU = BÇ(AÈA') = (BÇA) È (BÇA') volgens resp stelling 4, 7 en 15. De doorsnede is leeg want (BÇA) Ç (BÇA') = (BÇB) Ç (AÇA') = BÇÆ = Æ volgens stelling 11, 13 en 2. Dus P(B) = P(BÇA) + P(BÇA') = P(B|A).P(A) + P(B|A').P(A'). Een, maatschappelijk gezien, buitengewoon belangrijk gevolg van deze stelling is de notie dat de uitslag positief/negatief van een onderzoek niet voor honderd procent te vertrouwen is. Contraexpertise is absoluut aan te bevelen! Met name wat men noemt vals-positieve uitkomsten kunnen tot grote, onnodige ongerustheid leiden. Een onderzoek (bij een grote groep mensen) wordt altijd uitgevoerd met een test die niet voor 100% betrouwbaar (correct) is. Een voorbeeld. Aangenomen dat 1% van de mensen feitelijk besmet (ziek, gedrogeerd) is (B=0,01 is de fractie die besmet is) en aangenomen dat de test (het bloedonderzoek, de Röntgenfoto) in 99% van de gevallen een correcte uitslag geeft (C=0,99 is de fractie correcte uitslagen) dan is de geloofwaardigheid van een positieve uitslag (positief wil zeggen: "Helaas, U bent ziek/besmet/gedrogeerd"): G = P(besmet | test pos) = BC / [BC + (1-B)(1-C)] = 0,5 Met andere woorden: in dit geval hoeft iemand die een positieve uitslag krijgt ("Helaas, …") nog niet te wanhopen, want er is maar liefst 50% kans (hoop) op een gunstiger uitkomst, een vals-positieve uitkomst. Gauw contra-expertise aanvragen dus! De volgende tabel illustreert de verwachting volgens Bayes bij het bovenstaand voorbeeld (N=10 000, B=0.01, C=0.99). test correct
test fout
totaal aantal
besmet
B.C.N=99
B(1-C).N=1
B.N=100
niet besmet
(1-B).C.N=9801
(1-B)(1-C).N=99
(1-B).N=9900
totaal aantal
C.N=9900
(1-C).N=100
N=10 000
In het simulatieprogramma BAYES is dit uitgeprobeerd.
KANSEN
- 48 -
Wetenschappelijk? Soms halen wiskundige of statistische overwegingen de krant. Na het typen van de vorige regels las ik (in de NRC van 11 oktober 2003) het volgende in een ingezonden brief.
UMTS onderzoek Volgens een bericht over een onderzoek naar de mogelijke effecten van straling van mobiele telefoon zendmasten, zou er een 'statistisch significant' (d.w.z. met 95% betrouwbaarheid) verband zijn tussen de blootstelling aan straling van UMTS-antennes en het optreden van lichamelijk (on)welbevinden. Deze conclusie is fout. De onderzoekers zijn in een bekende (en beruchte) statistische valkuil getrapt. In werkelijkheid is het verreweg het waarschijnlijkst dat de uitkomst te wijten is aan puur toeval, en is er dus geen reden om UMTS-masten nu als 'verdacht' te bestempelen. Ik heb serieuze twijfels over de manier waarop het onderzoek is uitgevoerd. Maar aangenomen dat het onderzoek inderdaad een betrouwbaarheid heeft van 95%, is er dus 95% kans dat, als er in werkelijkheid geen verband is, er ten onrechte geconcludeerd wordt dat dat er wel is. Er is dan sprake van een vals-positieve uitkomst. De onderzoekers verbinden hieraan echter de verkeerde conclusie dat er dus met 95% zekerheid is vastgesteld dat blootstelling aan UMTS-straling waarneembare effecten oplevert. Dat dit fout is kan als volgt worden aangetoond. Ik schat de kans dat er inderdaad een effect is op (bijna) nul. Voorafgaand aan het onderzoek van TNO is er immers geen verband aangetoond. Ook geeft het TNO onderzoek zelf aan dat soortgelijke straling van GSM-antennes blijkbaar geen waarneembare effecten heeft. Laten we aannemen dat er bij voorbaat 1% kans is dat er inderdaad een effect is. De kans op een valspositieve uitkomst is dus 0.99´0.05 en is dus ongeveer 5 maal groter dan de kans van 0.01´0.95 dat het onderzoek de juiste uitkomst geeft. Dus is er maar liefst 84% kans dat de verkeerde conclusie getrokken is. Totaal anders dan de 5% onzekerheid gesuggereerd door de onderzoekers. Het bovenstaande ligt ten grondslag aan de formule van Bayes. De Tweede Kamer dient te overwegen of het nuttig is om nog meer belastinggeld te steken in een onderzoek waarvan het vrijwel onvermijdelijk is dat de uitkomst voornamelijk op toeval berust. Prof. Dr. Ir. B.J. van Wees, Groningen Ik maak hier een tabelletje bij. Dat komt voort uit een dwangmatige behoefte om uit te leggen.
conclusie:
effect geen effect totaal
werkelijk effect: ja nee 0.01´0.95 0.95´0.99 0.01´0.05 0.05´0.99 0.01 0.99
De kans op een verkeerde conclusie =
KANSEN
totaal 0.95 0.05 1
0.05 ´ 0.99 = 0.84 0.05 ´ 0.99 + 0.01 ´ 0.95
- 49 -
7
Extra opgaven over simulaties
23.
Laat zien dat in een klas van 30 leerlingen de kans dat er twee gelijk jarig zijn minstens 70% is. Doe dat op twee manieren: (a) Via 200 simulaties in een programma JAARSIM waarbij je 30 getallen zonder terugleg trekt uit U={1,2,3,…,30}. Sorteer deze 30 getallen (ASort) en onderzoek of de verschilrij (DList) nullen bevat. Zo ja, dan 30 nieuwe getallen genereren. Deze controle op het niet-terugleggen kost tijd waardoor je een beetje geduld moet hebben. Tweehonderd simulaties echter is wel het minste wat je kunt doen om duidelijkheid te krijgen. (b) De formule voor het verjaardagenprobleem in een klas van 30 leerlingen is:
1-
365 364 363 336 × × × ... × 365 365 365 365
Bedenk een procedure om de TI-83 dit te laten uitrekenen. Dat is niet zo gemakkelijk als je zou denken, want uitdrukkingen met een faculteit erin geven al gauw een OVERFLOW foutmelding. Misschien is het handig, om de MATH functie nPr te gebruiken. De functie N nPr K staat voor het aantal permutaties van rijen met K elementen, gekozen uit een totaal aantal van N elementen. Het is een soort machtsverheffen met teruglopende factoren: 10 nPr 3 geeft niet 10´10´10 maar 10´9´8=720. Neem, om dit toe te lichten, alle getallen van drie cijfers van 000 t/m 999. Hoeveel van die 1000 getallen bestaan uit drie verschillende cijfers? Het antwoord op die vraag is 10´9´8 = 10 nPr 3 = 720. Zet een formule in Y1 voor alle mogelijke klasgroottes X£40 in het verjaardagenprobleem en lees de tabel af die hoort bij Y 1. 24.
Het dronken beeldschermpunt. Programma RANDWALK. Startend links-midden op positie (30,5) beweegt pixel (X,Y) met stapgrootte 1 in principe naar rechts. Echter: het pixel is lichtelijk aangeschoten. De kansen op een stapje omhoog,omlaag en naar links hebben een (random) kans 0.2; het stapje naar rechts heeft een tweemaal zo grote kans (0.4 dus). Bijvoorbeeld. Als rand>0.6 doe je X+1® X; als rand<0.2 dan Y+1®Y; als rand tussen 0.2 en 0.4 dan Y-1®Y; als rand tussen 0.4 en 0.6 komt er X-1® X. Schrijf hiervoor een programmaatje. Pas op: na Pixl-On(X,Y) wordt de Y op nul gezet (de X niet, vreemd genoeg); dat is een fout in TIbasic die ook speelt bij de Pt-On en line functies. Je moet dus eerst Y even opzij zetten met Y®Z en na Pixl-On deze coördinaat weer oproepen met Z ®Y.
25.
Schrijf een simulatieprogramma BIN4KANT dat de gemiddelde afstand berekent tussen twee punten (A,B) en (C,D) binnen een eenheidsvierkant. Een eenheidsvierkant heeft de hoekpunten (0,0), (1,0), (1,1) en (0,1). Je kunt dus voor de coördinaten simpelweg de instructie rand nemen. De afstandsformule is volgens Pythagoras:
(C - A) 2 + ( D - B ) 2 26.
Dezelfde vraag voor twee punten binnen een eenheidscirkel. Let op 2 2 2 2 de voorwaarden: A + B < 1 en C + D < 1. Het programma heet BINCIRKL.
KANSEN
- 50 -
2
27.
De algemene kwadratische vergelijking ax + bx + c = 0 heeft, zoals je weet, geen 2 oplossing als de "discriminant" b - 4ac < 0 is. Je kunt dit ook anders formuleren: 2 De parabool y = ax + bx + c heeft geen snijpunt of raakpunt met de 2 x-as, als b - 4ac < 0. Zo te zien is de kans op het al of niet snijden van de x-as ongeveer fifty-fifty.(Het geval met maar één oplossing, het raken, is zo zeldzaam dat we het gevoeglijk kunnen vergeten). Maar …. Is dit weer een voorbeeld van gezichtsbedrog, zoals als we al veel vaker hebben meegemaakt of kunnen we dit keer op het gezonde verstand vertrouwen? We noemen het programma WORTELS2. Je kunt a, b en c random nemen tussen bijvoorbeeld -1 en 1 (of -10^10 en +10^10). Maar verstandiger is het, te beginnen met a=1, dus eerst te 2 kijken naar de wortels van x + bx + c = 0.
28.
a) Simuleer het werpen met een munt. Het verschil tussen het werkelijke aantal keren kop en het verwachte aantal kop neemt langzaam toe. Het verschil tussen de werkelijke fractie kop en de verwachte fractie kop (0.5) echter neem snel af naar nul. Laat dit zien in een grafiekje in het Draw scherm. Neem Ymin=-30 en Ymax=30, Xmin=0 en Xmax=94. Hierin is X het worpnummer en Y = (aantal kop - verwachte aantal kop). Zorg ervoor dat je de serie van 94 worpen telkens kunt herhalen, want 94 is te weinig om iets te kunnen zien van deze "stochastische wandeling". Onderaan het scherm is nog plek voor het verschil van de fracties. Daar zou je tegelijkertijd het fractieverschil -30+60´ fractieverschil kunnen weergeven. Dat zou naar nul moeten gaan. Je zult zien dat hoe groter het aantal worpen is, des te langer het duurt voor de eerstgenoemde stochastische wandeling terugkeert op nulniveau. Daar hebben we weer een voorbeeld van de grilligheid van de statistiek. Onze intuïtie heeft vaak moeite met zulke verschijnselen. Het programma hebben we STOWAND genoemd. b) Het duurt veel langer dan je zou denken voordat de gevonden fractie (de "zweetkans") een beetje gaat lijken op de theoretische kans (de "weetkans"). Hoeveel worpen (n) met een geldstuk heb je minstens nodig om de de-decimaal (d=1 is de eerste decimaal) met 95% waarschijnlijkheid goed te hebben? Geef een formule die n in d uitdrukt; bedenk dat in dit geval de standaarddeviatie s gelijk is aan Ö(n.p.q) dus Ö(n.½ .½) = ½Ön. en dat je volgens een vuistregel een afwijking 2s moet nemen.
KANSEN
- 51 -
ANTWOORDEN 1.
0.50
0.33 (vermoedelijk 1/3)
nee(1/3¹0)
2.
37 mensen hadden geen huisdier
ja (zie opg 17)
3. Met de kansen 0.5 en 0.5: ClrHome:ClrAllLists rand(400)üL1:rand(400)üL2 L1*L2üL3 1-Var Stats L1:ËüA 1-Var Stats L2:ËüB 1-Var Stats L3:ËüC Output(1,1,"Ë1=") Output(1,4,round(A,5) Output(2,1,"Ë2="):Output(2,4,round(B,5):Output(3,1,"Ë3=") Output(3,8,round(C,5):Output(4,1,"Ë1*Ë2="):Output(4,8,round(AB,5)) Met de kansen 0.5 en 0.8: ClrHome:ClrAllLists rand(500)üL1:rand(500)üL2 sum(L1<.5)üA:sum(L2<.8)üB sum(L1<.5 and L2<.8)üC round(A/500,5)üA round(B/500,5)üB round(C/500,5)üC Output(1,1,"Ë1=") Output(1,4,A):Output(2,1,"Ë2="):Output(2,4,B) Output(3,1,"x•ENx‚="):Output(3,8,C) Output(4,1,"Ë1*Ë2="):Output(4,8,round(AB,3)
Theorie
4.
5.
6. 7.
Meting
P( 0,495n £ X £ 0,505n | m=0,5n, s=0,5Ön ) = 0,95 oftewel: 0=normalcdf(0.495N,0.505N,0.5N,.5Ö(N))-0.95) levert N =38414 worpen
æ n ö 1 99 98 × × × ... = n. 1 = 0,01n ; hier is de kans dus 0,37. De kans is çç ÷÷ × 100 99 98 100 è1ø Alle personen hebben gelijke kans! De laatste bijvoorbeeld heeft de prijs, als al zijn 9 ×8×7×6×5×4×3×2×1 = 1 voorgangers gemist hebben: p = 10 ; bewonder het elegante 9 8 7 6 5 4 3 2 10 wegdelen van gelijke factoren in de tellers en noemers.
8.
KANSEN
- 52 -
b) 0=N*sum(seq(1/X,X,1,N))-2000 geeft N=316
9.
a) 519
10.
Volgens model II is de kans dat het andere kind ook een meisje is 1/3.
11.
a) 1+ 3+...+(2n-1) + (2n+1) = (n+1) ?? n + 2n + 1 = n + 2n + 1 klopt b) de halve rechthoek van 7 bij 8 is gevuld
12.
Na 3000 simulaties vond ik in L1 de lijst met gemiddelde 4.6667 volgens de theorie: met verwachting 4.6763 Haast te mooi, maar desalniettemin waar.
13.
Zonder terugleg: 2/12, 4/12 en 6/12
14.
1/3 n(n -1) + (n+1)n
15.
Voor n£6: recursieformule: Voor n³7 is de formule:
16.
Directe formule voor N£6: P(N=n)=1/6.(1+1/6) =7 /6 Een directe formule is voor N³7 niet gemakkelijk te geven. Wel is duidelijk dat P(N=n) voor grote waarden van n naar een limiet convergeert.
17.
De gemiddelde verwachting met een dobbelsteen is 3,5 = 7 2 ogen.
2
2
2
{0, 201, 390, 607, 783, 1019} {0, 200, 400, 600, 800, 1000}
Met terugleg: 1/16, 3/16, 5/16 en 7/16 2
= 1/3 (n+1)((n+1) -1) ?? = 1/3 (n+1)n(n+2) links en rechts maal 3/n geeft 2 n -1 + 3n+3 = (n+1)(n+2) ?? 2 2 n + 3n + 2 = n + 3n + 2 klopt P(N=n)=1/6.[1+P(N=1)+P(N=2)+…+P(N=n-1)] P(N=n)=1/6.[P(N=n-6)+P(N=n-5)+…+P(N=n-1)] n-1
De kans is op den duur dus 18.
2
1 7
n-1 n
= 2 7 » 0,2857
2
Er zijn 9 van de 24 gunstig. De kans is dus 9/24 of 3/8. De volgende formule geeft de juiste kans voor alle waarden van n: 1 1 1 1 1 - + - + ...( -1) n (in limiet staat hier 1/e). 2! 3! 4! 5! n! Dit kun je uitrekenen op de TI-83 met sum(seq((-1)^K/(K!),K,2,X))ü Y1
19.
20.
P(a<0.5)=0.75
21.
P(a<0.2)=0.64
KANSEN
- 53 -
22.
a) P(a<0.5)=2/3 b) De gearceerde oppervlakte verhoudt zich tot de totale
a
2
6 - (3 - 2a )(2 - 2a ) 5a - 2a = = 5 a - 2 a2 , 3 3 6 3 met a£1 (want P kan nooit verder dan 1 vanaf de zijden
a
oppervlakte als
liggen). De afgeleide hiervan is 53 - 43 a ; zodat de verwachte gemiddelde afstand als volgt te berekenen is: 1
a = ò a × ( 53 - 34 a ) da = 56 a 2 - 94 a 3 0
]
1 0
7 = 18 .
0üS:0üT Lbl 0 3randüX:2randüY min({X,Y,3-X,2-Y})üD T+1üT S+DüS Disp S/T Goto 0 23.
JAARSIM ClrHome 0üS For(X,1,200) randInt(1,365,30)üL1 SortA(L1) ¾List(L1)üL2 (sum(L2=0)ù1)+SüS Output(1,1,X) Output(2,1,S/X) End JAARBER "1-(365 nPr X)/(365^X)"üY1 0üTblStart:1ü¾Tbl DispTable
24.
RANDWALK ClrHome:ClrDraw:AxesOff FnOff :PlotsOff 30üX:5üY Lbl 0 YüZ "Y COORD TIJDELIJK" Pxl-On(X,Y) "OPSLAAN" ZüY "Y WEER OPHALEN" randüR If R<.2:X+1üX If Rù.2 and R<.4:X-1üX If Rù.4 and R<.6:Y-1üY If Rù.6:Y+1üY Goto 0
KANSEN
- 54 -
2
2-2a 3-2a 3
P
25.
De gemiddelde afstand is ongeveer 0.52 BIN4KANT ClrHome:0üS For(X,1,200) randüA randüB randüC randüD ð((C-A)Ü+(D-B)Ü)üE S+EüS Output(1,1,S/X) End
26.
De gemiddelde afstand is ongeveer 0.91 BINCIRKL ClrHome 2üA:2üB:2üC:2üD:0üS For(X,1,400) Repeat AÜ+BÜ<1 2rand-1üA 2rand-1üB End Repeat CÜ+DÜ<1 2rand-1üC 2rand-1üD End ð((C-A)Ü+(D-B)Ü)üE S+EüS Output(1,1,S/X) End
KANSEN
- 55 -
27.
De kans op een oplossing lijkt ongeveer 62,6% ClrHome:0üT For(X,1,500) 2rand-1üA:2rand-1üB:2rand-1üC If BÜ-4AC>0:T+1üT Output(1,1,T/X) End Het geval a = 1: Het gaat nu om de kans P op een niet-negatieve discriminant 2 b – 4c ³ 0 met b en c random tussen de grenswaarden –G en +G: G*(2rand-1) dus. Bekijk het vierkant met zijde 2G en daarin drie grafieken van 2
c = ¼ b , voor verschillende waarden van G. De c-as vertikaal, de bas horizontaal. Voor G = 4 gaat zo'n grafiek door het hoekpunt (G, G). Voor G < 4 is de kans op oplossingen te vinden door de gearceerde G
oppervlakte
ò ( 14 b
2
- - G ) db =
1
6G
3
+ 2G 2 te delen door de
-G 2
oppervlakte 4G van het vierkant. Dat levert voor G < 4 dus P = 1 G + 1 ; als G=1 genomen wordt beschrijf je voorgaande 24 2 simulatie en komt er
13 » 0,5417 uit. 24
Voor G ³ 4 is de gearceerde oppervlakte gelijk aan 2 G
ò ( 14 b
2
2 - - G ) db + (G - 2 G ) × 4G ; delen door 4G geeft P = 1-
-2 G
2 3 G
.
Even uitproberen met 16(2rand-1) voor B en C. Voor deze waarde van G (=16) zou de kans P = 0,8333 moeten zijn. Ik vond 0.8334225721 na 69000 simulaties. Niet gek dus. Heel verrassend is dat P naar 1 gaat voor grote waarden van G, hetgeen betekent dat de 2 kans dat een parabool y = x + bx + c de x-as snijdt vrijwel 100% is, voor b,c Î Â. Wie dat begrijpt, mag het zeggen. Het heeft natuurlijk iets te maken met het domein van b en c: voor grote waarden van G telt het witte binnengebied van de parabool nauwelijks meer mee. 2
2
Het algemene geval y = ax + bx + c met b – 4ac ³ 0 en a,b,c Î Â schijnt een ander antwoord, niet minder raadselachtig, op te leveren. Na een groot aantal simulaties met verschillende G-waarden kreeg ik het vermoeden dat P, onafhankelijk van het domein van de randomwaarden, ongeveer 0,626 is. Maar dat zou niet moeten kunnen, omdat de ontaarde waarde a = 0 roet in het eten zou moeten gooien. Zou moeten ... En waarom moet ik nu ineens aan de Bone Man denken?
KANSEN
- 56 -
28.
a) STOWAND Lbl 0 For(X,1,94) S+randInt(0,1)üS:X+94TüU round(.5U,0)üV:S-VüY:Pt-On(X,Y) If fPart(U/10)=0:Text(0,32,U) abs(Y/U)üF:ú30+60FüG:Pt-On(X,G) End T+1üT:Goto 0 –d –d 2d b) 2s/ m = 10 geeft 2Ön/n= 10 dus n = 4.10 . Dat wil zeggen, voor een goede (betrouwbare) schatting van de derde decimaal moet je minstens 4 miljoen simulaties uitvoeren!
KANSEN
- 57 -
