Card Problem from Gardner

Question

The following problem is from aha! Insight by Martin Gardner.

On a draw of $n$ cards from a deck without replacement, what is the probability that seven of the $n$ cards all have the same suit?

It has been a while since I have done a probability computation so I am here to have my solution verified.

Let $X_n$ be the random event of drawing at least seven cards of the same suit upon drawing $n$ cards.

$$P(X_n) = \left\{ \begin{array}{ll} 0 & : n < 7\\ \dfrac{\binom{4}{1}\binom{13}{7}\binom{39}{n-7}}{\binom{52}{n}} & : 7 \le n < 25\\ 1 & : 25 \le n \end{array} \right. $$

EDIT: So this counts the number of ways to draw exactly seven cards with $n$ draws. Robert's post gives the desired answer, but I am wondering if there is a simpler expression for the probability given solely in terms of $n$.

EDIT2: Started bounty looking for a solution in terms of $n$. I would also accept an explanation as to why Robert's solution is as simple as it gets.

Answer 1

Let's consider a simpler problem with only two suits, say spades and hearts. Lets suppose there are $S$ cards of each suit in the deck, that you will choose $n$ cards, and that you want to know the probability of getting $k$ or more cards of at least one suit.

The probability of getting exactly $k$ spades is $$ \frac{\binom{S}{k}\binom{S}{n-k}}{\binom{2S}{n}}. $$ The probability of getting $k$ or more spades is $$ \sum_{j=k}^n\frac{\binom{S}{j}\binom{S}{n-j}}{\binom{2S}{n}}. $$ The probability of getting $k$ or more cards of at least one suit can be obtained using the principle of inclusion-exclusion by adding the probability of getting $k$ or more spades to the probability of getting $k$ or more hearts, and then subtracting the probability of getting $k$ or more of both suits. This computation gives $$ 2\sum_{j=k}^n\frac{\binom{S}{j}\binom{S}{n-j}}{\binom{2S}{n}}-\sum_{j=k}^{n-k}\frac{\binom{S}{j}\binom{S}{n-j}}{\binom{2S}{n}}. $$ There is no simple, closed form for these sums. They can be expressed in terms of the cumulative distribution function of the hypergeometric distribution. The linked Wikipedia article gives the formula for the cumulative distribution function in terms of a ${}_3F_2$ generalized hypergeometric function, but this is no help, since the generalized hypergeometric function is defined to be precisely the kind of sum of products that I presume you're trying to get away from.

There is an algorithm for determining when sums such as these have a closed form, under suitable definition of "closed form". See the Wikipedia articles on Gosper's algorithm and Wilf–Zeilberger pairs and the references in those articles for information about this. I am not in a position at the moment to tell you what this theory has to say about the sums above, but, suffice it to say, there is no closed form.

When you increase the number of suits to four, things get more complicated. The principle of inclusion-exclusion yields the formula $$ \begin{aligned} &4\sum_{h=k}^n\frac{\binom{S}{h}\binom{3S}{n-h}}{\binom{4S}{n}}-6\sum_{h=k}^n\sum_{i=k}^{n-h}\frac{\binom{S}{h}\binom{S}{i}\binom{2S}{n-h-i}}{\binom{4S}{n}}+4\sum_{h=k}^n\sum_{i=k}^{n-h}\sum_{j=k}^{n-h-i}\frac{\binom{S}{h}\binom{S}{i}\binom{S}{j}\binom{S}{n-h-i-j}}{\binom{4S}{n}}\\ &\qquad-\sum_{h=k}^n\sum_{i=k}^{n-h}\sum_{j=k}^{n-h-i-k}\frac{\binom{S}{h}\binom{S}{i}\binom{S}{j}\binom{S}{n-h-i-j}}{\binom{4S}{n}} \end{aligned} $$ for the probability of $k$ or more cards of at least one suit. Only the first of the sums in this formula can be expressed in terms of generalized hypergeometric functions. I am not sure whether the Wilf-Zeilberger method can be extended to multiple sums, but, given the difficulties that present themselves even in the two-suit case, it is almost certain that there is no closed form.

Answer 2

For each four-tuple $(c,d,h,s)$ of nonnegative integers with $c+d+h+s = n$ and $\max(c,d,h,s) \ge 6$, the probability of getting $c$ clubs, $d$ diamonds, $h$ hearts and $s$ spades is $$ p(c,d,h,s) = \dfrac{{13 \choose c}{13 \choose d} {13 \choose h} {13 \choose s}}{{52 \choose n}}$$ You then have to add this over all $(c,d,h,s)$. For example, if $n= 15$ the possible $(c,d,h,s)$ are $$[2, 3, 3, 7], [2, 2, 4, 7], [1, 3, 4, 7], [0, 4, 4, 7], [1, 2, 5, 7], [0, 3, 5, 7], [1, 1, 6, 7], [0, 2, 6, 7], [0, 1, 7, 7], [2, 2, 3, 8], [1, 3, 3, 8], [1, 2, 4, 8], [0, 3, 4, 8], [1, 1, 5, 8], [0, 2, 5, 8], [0, 1, 6, 8], [0, 0, 7, 8], [2, 2, 2, 9], [1, 2, 3, 9], [0, 3, 3, 9], [1, 1, 4, 9], [0, 2, 4, 9], [0, 1, 5, 9], [0, 0, 6, 9], [1, 2, 2, 10], [1, 1, 3, 10], [0, 2, 3, 10], [0, 1, 4, 10], [0, 0, 5, 10], [1, 1, 2, 11], [0, 2, 2, 11], [0, 1, 3, 11], [0, 0, 4, 11], [1, 1, 1, 12], [0, 1, 2, 12], [0, 0, 3, 12], [0, 1, 1, 13], [0, 0, 2, 13], [0, 0, 1, 14], [0, 0, 0, 15] $$ and their rearrangements.

In Maple:

with(combinat):
for n from 7 to 24 do
  P:= map(`+`,select(t -> nops(t) = 4 and max(t) >= 8, partition(n+4)),
          [-1,-1,-1,-1]);
  prob:= add(mul(binomial(13,p[i]),i=1..4)
          *nops(permute(p))/binomial(52,n),p=P);
  printf("The probability for n = %d is %a, or approximately %4.5f\n",
         n, prob, evalf(prob));
od:

The probability for n = 7 is 33/643195, or approximately 0.00005

The probability for n = 8 is 1166/3215975, or approximately 0.00036

The probability for n = 9 is 272/189175, or approximately 0.00144

The probability for n = 10 is 116534/27657385, or approximately 0.00421

The probability for n = 11 is 40148/3951055, or approximately 0.01016

The probability for n = 12 is 40654/1905803, or approximately 0.02133

The probability for n = 13 is 44677/1108025, or approximately 0.04032

The probability for n = 14 is 397745888/5669763925, or approximately 0.07015

The probability for n = 15 is 491355875/4309020583, or approximately 0.11403

The probability for n = 16 is 139443061606/797168807855, or approximately 0.17492

The probability for n = 17 is 11953182762/46892282815, or approximately 0.25491

The probability for n = 18 is 11865928058/33494487725, or approximately 0.35427

The probability for n = 19 is 14099080786/29968752175, or approximately 0.47046

The probability for n = 20 is 3579684638/5993750435, or approximately 0.59724

The probability for n = 21 is 60779104513/83912506090, or approximately 0.72432

The probability for n = 22 is 1090469390477/1300643844395, or approximately 0.83841

The probability for n = 23 is 261904914329/282748661825, or approximately 0.92628

The probability for n = 24 is 8032961212957/8199711192925, or approximately 0.97966

Answer 3

Here is a Mathematica code that returns Robert's computations above:

Table[f[list_] := Binomial[13, list[[1]]] Binomial[13, list[[2]]] Binomial[13, list[[3]]] Binomial[13, list[[4]]]/Binomial[52, n]; Total[Map[f, Select[Compositions[n, 4], Max[#] >= 7 &]]], {n, 7, 25}]