Vtyjkyuk

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.

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 $binom13k$ ways to select $k$ particular conditions, so by generalized inclusion-exclusion there are

$$
sum_k=1^13(-1)^k-1binom k1binom13k(52-3k)!(4!)^k
$$

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

$$
frac26678312213573602672321879823867555861136785646174010387823491412456845703125approx0.00234682
$$

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

$$
frac13cdot49!cdot4!52!=frac13cdot2452cdot51cdot50=frac1425approx0.00235294;.
$$

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

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

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.