What is the probability of dealing all $52$ cards in a standard well shuffled deck getting a single quad?

Question

I would like to know if you take a well shuffled deck of $52$ cards and then deal out all of them one at a time without replacement, what is the probability that they will be dealt such that there is exactly one quad dealt in order such as $4$ consecutive kings? That run of $4$ has to be the only run of that length, there cannot be any others. For example, 5, ... ,K,K,K,K,7,J,J,J,J is a "loser". Other than computer simulation I am not sure how to solve this.

Answer 1

This can be done by inclusion-exclusion. We have $13$ conditions: the existence of a run of one of the $13$ card values. The number of configurations in which at least $k$ particular of these conditions are fulfilled is $(52-3k)!(4!)^k$, since we can consider the runs of $4$ as atoms in shuffling the cards, thus reducing the number of items to be shuffled by $3$ per run, and then we have $4!$ permutations within each run.

There are $\binom{13}k$ ways to select $k$ particular conditions, so by generalized inclusion-exclusion there are

$$ \sum_{k=1}^{13}(-1)^{k-1}\binom k1\binom{13}k(52-3k)!(4!)^k $$

configurations with exactly one run. Dividing by the total number $52!$ of configurations yields the probability

$$ \frac{2667831221357360267232187982386755586}{1136785646174010387823491412456845703125}\approx0.00234682 $$

of getting exactly one run. Since multiple runs are unlikely, the result is already quite well approximated by the first term,

$$ \frac{13\cdot49!\cdot4!}{52!}=\frac{13\cdot24}{52\cdot51\cdot50}=\frac1{425}\approx0.00235294\;. $$

Answer 2

Joriki provides the elegant solution ... but also shows the value of checking theoretical work via simulation. The discussion got me playing with simulating it in Mathematica, and thought I would post some code.

A pack of cards can be represented by the set:

cards = Table[Range[1, 13], 4] // Flatten 

Sampling without replacement, and dealing out a hand of 52 cards i.e. RandomSample[cards] ... and then counting all the cases where there are four identical cards in a row {x_,x_,x_,x_}, and doing it 10 million times:

4Row = ParallelTable[Count[Split[RandomSample[cards], SameQ], {x_,x_,x_,x_}], {10^7}];

All done. The empirical distribution is given by:

Tally[4Row]

{{0, 9976557}, {1, 23415}, {2, 28}}

In other words, out of 10 million drawings, there was:

• exactly 1 sequence of 4 identical cards 23,415 times (out of 10 million)

• exactly 2 sequences of 4 identical cards just 28 times (out of 10 million)

The above Monte Carlo simulation (10 million times) takes about 42 seconds on my Mac. I am sure there are faster ways of doing it, but that's another question :)

Answer 3

It seems that your shuffling algorithm is too slow. For the first card, you can pick a number n between 1 and 52. If it's 1, your card is where it's supposed to be. Otherwise, swap cards 1 and n. Then for card 2, pick a number between 2 and 52... While swapping seems to double the amount of movement, it means that the card to be chosen comes from a tight array so you don't have to search for it and this will make shuffling much faster.

For detecting 4 in a row, my idea is like punching holes in data cards so that if light shines through 4 cards, it's a match. The digital analog to this is setting a bit according to the rank of the card and ANDing four consecutive values together to see if one bit remains set.

Processed a billion simulations in 344 seconds with 1, 2, and 3 matches occurring with frequency 2347407, 3090, and 6 respectively.

program four
   implicit none
   real harvest(51)
   integer i, j
   integer, parameter :: N(51) = [(i,i=52,2,-1)]
   integer indices(51)
   integer matches(13)
   integer jmax
   integer cards(52)
   integer, parameter :: init(52) = [((ishft(1,j),j=0,12),i=1,4)]
   integer match
   integer(selected_int_kind(18)) t0,tf,rate

   matches = 0
   jmax = 1.0e9
   call random_seed()
   call system_clock(t0)
   do j = 1, jmax
      call random_number(harvest)
      indices = harvest*N
      cards = init
      do i = 1, 51
         if(indices(i) /= 0) cards([i,indices(i)+i]) = cards([indices(i)+i,i])
      end do
      match = 0
      do i = 1, 49
         if(iall(cards(i:i+3)) /= 0) match = match+1
      end do
      if(match > 0) matches(match) = matches(match)+1
   end do
   write(*,'(a,1x,i0)') 'Number of simulations = ',jmax
   write(*,'(a)') ' N Frequency'
   do i = 1, 13
      write(*,'(i2,1x,i0)') i , matches(i)
   end do
   call system_clock(tf)
   call system_clock(count_rate=rate)
   write(*,'(a,1x,g0)') 'Time in seconds = ', real(tf-t0)/rate
end program four

Answer 4

I wanted to post my algorithm for those that may want to attempt to simulate this in any general purpose language. These are the major steps . Some minor details are left out for simplicity here...

Step $1$: declare and initialize an array A of size $52$ such that A[$1$] = $1$, A[$2$] = $2$... A[$52$] = $52$.

Step 2: declare an initialize and array F[$52$] of flags to false so we know which cards are already drawn/picked. Change the appropriate flag to true when a card is picked (such as A[$22$] = true) if random number generator gives us $22$.

Step 3. Convert all $52$ cards to ranks only ($1$ to $13$)

Step 4: Count up runs of $4$. I "cheat" when I do this by converting the cards to a string such as "K2J7777A5..." then I just scan for AAAA, 2222, 3333... KKKK. If I detect exactly one set of quads it is a winning hand.

My first tweak is I only choose $51$ cards randomly and then just assign the one unchosen card as the last. That sped things up by about $10$%.

My next tweak is to try a faster way to choose all $52$ cards (such as simulate a shuffle). The language I am using has an array element delete function that will compact the remaining array so I can pick smaller and smaller random numbers that way but I don't yet know if it is faster.