Probability problem -- Rube Goldberg solution?

Question

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?

Answer 1

Don't even know how I came across this old question, but, interesting and seeing as Mr. W asked it, it piqued my interest.

The answer by David is a good use of simulation, and his "back-o-the-envelope" approximation using the probability functions of Mathematica over a discrete uniform is also kind of neat.

The other two answers, however are flat-out wrong: Both admit hand sets that are impossible, e.g., looking at the first entry of allPossibilities from Rojo's answer shows an entry of {{1, 1}, {1, 1}, {1, 1}}. I don't know about you, but "a standard 52 card deck" with six aces would get you shot in the old west...

In addition, counting winning hands via tuples/etc. and getting a ratio of winning to total hands is not correct (although by happenstance, it gets pretty close): take, e.g., the hand sets {{4, 4}, {3, 3}, {3, 3}} vs {{6, 5}, {4, 3}, {2, 1}}. In both, the "first" hand wins. But the second set is over 14X more likely to occur. So any tuple-counting/filtering approach must calculate the probabilities of each winning set and total them...

So, why not solve it directly and exactly?

(* bulid lists of possible ways to take from ranks of deck in two draws, pairs and singlets *)

two = Permutations[Prepend[ConstantArray[0, 12], 2]];
one = Permutations[Join[{1, 1}, ConstantArray[0, 11]]];
tt = Tr /@ (two*Range@13);
oo = Tr /@ (Range@13*# & /@ one);
ways = Sort@Join[Transpose[{oo, one}], Transpose[{tt, two}]];

(* helper functions - get p of a draw sum given current deck configuration,
   and p of draw sums less than given sum *)
pe[n_, db_] := With[{c = Select[ways, #[[1]] == n &]},
  Module[{z = #, r}, 
     If[FreeQ[r = db - z, _?Negative], {PDF[
        MultivariateHypergeometricDistribution[2, db], z], r}, {0, r}]] & /@ c[[All, 2]]]

ple[n_, db_] := With[{c = Select[ways, #[[1]] < n &]},
  Module[{z = #, r}, 
     If[FreeQ[r = db - z, _?Negative], {PDF[
        MultivariateHypergeometricDistribution[2, db], z], r}, {0, r}]] & /@ c[[All, 2]]]


(* get probabilities of paths for given sum and remaining 2 players both having less *)
Monitor[res = Table[{n,
     peres = pe[n, ConstantArray[4, 13]];
     secondres = Map[ple[n, #] &, peres[[All, 2]]];
     step3 = Transpose[{secondres[[All, All, 1]]*peres[[All, 1]], secondres[[All, All, 2]]}];
     step4 = Map[With[{c = #}, Map[ple[n, #] &, c]] &, step3[[All, 2]]];
     Total[Total /@ Map[Total, step3[[All, 1]]*step4[[All, All, All, 1]], {2}]]
     }, {n, 2, 26}], n];

(* final probability of a player winning *)

res[[All, 2]] // Total
% // N

(*
4710593/15268890
0.308509
*)

And some graphical results (a priori probability of holding some sum, probability of sum given you won, and probability of winning given you are holding some sum):

Answer 2

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

Answer 3

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.

Answer 4

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.

Answer 5

f[u_] := And[#1 > #2, #1 > #3] & @@ (Plus @@@ Partition[Mod[u, 13, 1], 2])
fun[n_] := With[{r = Table[RandomSample[Range[52], 6], {n}]}, 
           N[Length@Pick[r, f /@ r, True]/n]]

Not particularly efficient but:

fun[10000000]

yielded: 0.308507