/* SAS program to accompany the blog post "Head-tail versus head-head: A counterintuitive property of coin tosses" by Rick Wicklin, published 13APR2016 http://blogs.sas.com/content/iml/2016/04/09/counterintuitive-property-of-coin-tosses */ /* Story from Wired magaxine: http://www.wired.com/2016/03/mathematicians-discovered-prime-conspiracy/ "A lecture at Stanford by the mathematician Tadashi Tokieda, of the University of Cambridge, in which he mentioned a counterintuitive property of coin-tossing: If Alice tosses a coin until she sees a head followed by a tail, and Bob tosses a coin until he sees two heads in a row, then on average, Alice will require four tosses while Bob will require six tosses (try this at home!), even though head-tail and head-head have an equal chance of appearing after two coin tosses." */ data CoinFlip(keep= player toss); call streaminit(1234567); label Player = "Stopping Condition"; do i = 1 to 100000; do player = "HT", "HH"; /* first Alice tosses; then Bob */ first = rand("Bernoulli", 0.5); /* 0 is Tails; 1 is Heads */ toss = 1; done = 0; do until (done); second = rand("Bernoulli", 0.5); toss = toss + 1; /* how many tosses until done? */ if player = "HT" then /* if Alice, ... */ done = (first=1 & second=0); /* stop when HT appears */ else /* if Bob, ... */ done = (first=1 & second=1); /* stop when HH appears */ first = second; /* remember for next toss */ end; output; end; end; run; proc means data=CoinFlip MAXDEC=1 mean median Q3 P90 max; class player; var toss; run; /* comparative histograms */ ods graphics / height=500px width=600px; title "Distribution of Coin Tosses until HT or HH Observed"; title2 "100,000 Games Simulated"; proc sgpanel data=CoinFlip; panelby Player / rows=2 layout=rowlattice rowheaderpos=left; histogram Toss / binwidth=1; colaxis values=(2 to 30); rowaxis grid; run; /* comparative ECDFs */ proc univariate data=CoinFlip; label Player = "Stopping Condition"; class player; cdfplot toss / overlay vref=50 75 90; ods select cdfplot; run; /* simulate 100,000 games. Keep track of how many Alice and Bob win */ data Game(keep=winner); call streaminit(1234567); do i = 1 to 100000; first = rand("Bernoulli", 0.5); /* 0 is Tails; 1 is Heads */ done = 0; do until (done); second = rand("Bernoulli", 0.5); if (first=1 & second=0) then do; done = 1; winner = "HT"; /* stop when HT appears */ end; if (first=1 & second=1) then do; done = 1; winner = "HH"; /* stop when HH appears */ end; first = second; /* remember for next toss */ end; output; end; run; proc freq data=Game; tables winner / nocum binomial(level='HT'); run;