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):
