# Probability problem — Rube Goldberg solution?

A user posted this question on StackOverflow which was closed as off topic:

3 people are playing a game with a standard 52 card deck. Each player is given 2 cards each, possible cards and their values being $\mathrm A = 1$, $2=2$, $\ldots$, $\mathrm J=11$, $\mathrm Q=12$, $\mathrm K=13$. A player wins if the sum of his cards is greater than both of the sums of the other player's cards. What is the winning probability for each player?

I am curious to see this problem receive the kind of solution that David gave here. Allowing Mathematica to do as much of the thinking for you as possible, how can this be solved?

-

Copying my clock post? Impossible! Anyway,

# Chapter one: using brute force.

Before you complain: I'm a physicist, this is how we do mathematics: by experiment. Ha!

We need cards! Inconveniently, Mathematica currently lacks built-in support for Kings. We therefore have to use a workaround: let's call the named cards by their numbers. A is 1, J is 11 etc. The deck of cards consists of the range of 1-13 four times, giving the deck array

deck = ConstantArray[Range[13], 4] // Flatten

{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13}


Now let's setup the win/lose balance. While seemingly complicated in structure, this can be boiled down to a 2-tuple,

{win, lose} = {0, 0};


It is now time to ask the infinite monkeys to stop typing for a second and help us out play a couple of games, say ten million of them.

Do[picks = RandomSample[deck, 6];
If[Total@picks[[ ;; 2]] > Total@picks[[3 ;; 4]] \[And]
Total@picks[[ ;; 2]] > Total@picks[[5 ;;]],
++win,
++lose
];,
{10000000}
];


Run it! (Note: Monkeys don't seem to like being run in parallel.)

(Music entertains us while waiting)

Thanks guys, back to writing Shakespeare.

As for us, I've always been looking for an opportunity for making good use of BarChart3D.

Row@{"The probability of winning is ", 100 win/(win + lose) // N, " %."}
BarChart3D[{win, lose}]

The probability of winning is 30.8485 %.


# Chapter two: using bruter force.

If the previous solution isn't accurate enough, we can use Mathematica's statistics sledge hammer, which happily kills all high school probability problems. In the following I'm assuming that the deck of cards is infinitely large (uniform distribution), which is close enough for all practical purposes.*

dist = DiscreteUniformDistribution[{1, 13}];
Probability[a + b > c + d \[And] a + b > e + f, Thread[{a, b, c, d, e, f} \[Distributed] dist]]


$\displaystyle\frac{8800}{28561}~(=30.8112\;\%)$

That's close enough to the experimental result, at least to astronomy standards.

*: Straight out lie to cover up that I haven't thought of how to determine the actual distribution

-

Ok, I'll stick to "making MMA do the thinking"

The list of all the possibile hands a player can have

playerPossibilities = Tuples[Range[13], {2}];


The list of all the triplets of hands the 3 players could have been handed

allPossibilities = Tuples[playerPossibilities, {3}];


A function that returns True if a1 won

a1winsQ[{h1 : {_, _}, h2_, h3_}] := # > Total[h2] && # > Total[h3] &@Total[h1]


Count favourable cases versus total cases

Count[allPossibilities, _?a1winsQ]/Length[allPossibilities]

8800/28561

0.308112


OR

Perhaps you were looking for something more like this

In[46]:= cards = {"A", "2", "3", "4", "5", "6", "7", "8", "9", "10",
"J", "Q", "K"};
giveHands[n_: 1] := RandomChoice[cards, {n, 2}];

In[49]:= values = Thread[cards -> Range[13]]

Out[49]= {"A" -> 1, "2" -> 2, "3" -> 3, "4" -> 4, "5" -> 5, "6" -> 6,
"7" -> 7, "8" -> 8, "9" -> 9, "10" -> 10, "J" -> 11, "Q" -> 12,
"K" -> 13}

In[53]:= a1winsQ[{h1 : {_, _}, h2_, h3_}] :=
Unevaluated[# > Total[h2] && # > Total[h3] &@Total[h1]] /. values

play[] := a1winsQ[giveHands[3]];

In[65]:= res = 0;
Dynamic[N@res/i]
Do[res += Boole[play[]], {i, 1000000}]
res/1000000.

-
Tally[
Block[{me = #[[1]] + #[[2]]},
And[me > #[[3]] + #[[4]],
me > #[[5]] + #[[6]]]] & /@ (Permutations[
ConstantArray[Range[13], 4] // Flatten, {6}])]

(* ==> {{False, 3338816}, {True, 1487044}} *)

N[1487044/(3338816 + 1487044)]

(* ==> 0.308141 *)


edit a little better and with a little parallelization:

(Count[#, True]/Length[#]) &[
ParallelMap[
Block[{me = #[[1]] + #[[2]]},
And[me > #[[3]] + #[[4]], me > #[[5]] + #[[6]]]] &,
Permutations[ConstantArray[Range[13], 4] // Flatten, {6}]]]

(* ==> 28597/92805 *)

% // N

(* ==> 0.308141 *)


In both cases, we list all possible ways to deal out 6 of the cards (represented by 4 copies of 1–13), and for each way, determine whether or not player 1 wins, then count the number of ways player 1 wins and divide by the total number of ways of dealing out the 6 cards.

-