Friday, October 6, 2017

Card tricks brute forced

(Now going from lowest-level to highest-level coding)


When I was in middle school, we used a standard 36-card playing deck (from 6 onwards) to play a bunch of games. We also used this deck to do some "loves me, loves me not" type fortune telling where "loves me" was indicated by, essentially, an occurrence of two aces in a row somewhere in the deck.

Much later on I got curious about how likely that outcome was. A short back of the envelope calculation gives
P = [P(A1)*P(A2|A1)]*N 
where 
P(A1) = probability of "there is an ace at position 1 in the deck", 1/9
P(A2|A1) = probability of "there is an ace  at position 2 n the deck if there already is an ace at position 1", 3/35 (3 aces left, out of 35 cards possible)
N = 35 -- number of different positions in the deck where two consecutive aces can occur.

Hence

P = [(1/9)*(3/35)]*35 = 1/3


This looks too high, at first sight. 
So let's write a simple brute-force proof in Wolfram Language (ex-Mathematica):
cards = "6789TJQKA";
adeck = Characters[cards<>cards<>cards<>cards];
testdeck[deck_List] := Length[StringPosition[StringJoin@@ToString/@deck,"AA"]]>0;
(*Monte-Carlo*)
n=100000;For[c=0;i=0,i<n,i++,(thedeck=RandomSample[adeck];If[testdeck[thedeck],c++])]
approx=N[c/n]

It is more or less obvious what the code does: it just tests a large number n of permutations on a deck (using RandomSample[])and then testing whether the resulting deck has at least two aces in a row (searching for "AA" in a deck converted to a string) and counting the number of decks in which this is found.

You can run the code in Wolfram Cloud and make sure the answer is indeed quite close to 0.3
(In hindsight, is it so strange, really? A "loves me, loves me not" trick will only survive among kids if it has a significantly non-negligible chance of a favorable outcome. And I do believe this is how most fortune telling works.)


Naturally you can modify the code very simply to include other card tricks.
For example, you can set cards="23456789TJQKA" for a standard 52-card deck; additionally you can append <>"XX" to cards<>cards<>cards<>cards to include jokers (or <>"AA" if you want them to be wild and count as aces).
Or you can change the testing function to test for any card sequence (like "AK" or "6789"), or indeed modify the criterion to any degree of complexity.

Finally, if you want to differentiate the suits of cards, we will need a slightly more complicated code:

cards = "6789TJQKA"; suits = "SHDC";
adeck = Outer[StringJoin, cards//Characters, suits//Characters]//Flatten;
testdeck[deck_List] := Length[StringPosition[StringJoin@@(First[Characters[#]]&)/@deck,"AA"]]>0;
Isn't it just beautiful how you can combine Map (/@) and Apply (@@) along with Mathematica's pure functions to create tests on lists without unnecessary boilerplate like loops? (And we can verify that making all aces distinct does not change the probability of two of them following each other.)

Finally, I initially planned to do a rigorous (rather than Monte-Carlo) test using Permutations[adesk], but it turns out that Mathematica won't have lists with "more than machine integer" elements. Sad but would probably work (and chances are, will be much faster than Monte-Carlo) with shorter decks.

No comments:

Post a Comment