How to find the probability of winning this game?

Question

Here's a very simple "game":

• Standard 52 card deck, shuffled.

• You deal cards face up, one at a time, while counting up towards 13, starting back at 1 once you reach 13.

• If the number you count out is the number of the card, you lose, reshuffle, and start over.

If you get through the whole deck, counting to 13 four times without matching up your counted number with the flipped card number, you win.

---

I programmed the game and ran it a few hundred million times. It seems like the probability of winning is about 1.65% (unless my game logic is flawed!).

How would I manually calculate the probability of winning?

Answer 1

This is tricky to solve because the events "the number you count out at step $i$ is the number on the card" are not independent.

Wikipedia, in its article on derangements, cites an interesting formula for the solution. It yields the answer

$$\frac{4610507544750288132457667562311567997623087869}{284025438982318025793544200005777916187500000000} \approx 0.016\,232\,727\,467\,194\,\ldots$$

---

A word is a sequence of letters from an alphabet. A derangement of a word is a permutation of its letters that differs from the original word in every place. The question concerns a 52-character word constructed from an alphabet of 13 distinct "letters" (the cards). The shuffled deck is the permutation and the game is won only when the permutation is a derangement.

The formula counts the number of derangements of a word having $r$ distinct letters whose counts are $n_1, n_2, \ldots, n_r.$ Let's call this $D(n_1,\ldots, n_r).$ The formula is

$$D(n_1,\ldots,n_r) = \pm \int_0^\infty \mathcal{L}_{n_1}(x)\mathcal{L}_{n_2}(x)\cdots\mathcal{L}_{n_r}(x)\,e^{-x}\mathrm{d}x$$

where the $\mathcal{L}_{n_i}$ are the Laguerre polynomials. In the question $r=13$ and all the $n_i=4,$ where

$$\mathcal{L}_4(x) = (x^4 - 16x^3 + 72x^2 - 96 x + 24)/4!.$$

Thus, "all" we need to compute is

$$D(4,4,\ldots, 4) = \frac{\pm 1}{(4!)^{13}}\int_0^\infty (x^4 - 16x^3 + 72x^2 - 96 x + 24)^{13} e^{-x}\mathrm{d}x.$$

That's actually straightforward, because after expanding the power in the integrand a polynomial $$\eqalign{p(x) &= x^{52} + p_{51}x^{51} + \cdots + p_1 x + p_0\\&= x^{52} - 208 x^{51} + 20904 x^{50} - 1352416 x^{49} + \cdots + 876488338465357824}$$ appears and the integral becomes a linear combination of Gamma functions (factorials, essentially) with value

$$D(4,4,\ldots, 4) = \frac{\pm1}{(4!)^{13}}(52! + p_{51}51! + \cdots + p_1 1! + p_0 0!).$$

It's easiest to use a computer to find the coefficients $p_i.$ You need to perform the addition to high precision due to the cancellation of large values that occurs (notice the alternating signs of the coefficients).

To obtain a probability, invoke the implicit assumption that all permutations are equally likely (a fair shuffle). The number of permutations of the word is

$$P(n_1,n_2,\ldots, n_r) = \binom{n_1+n_2+\cdots+n_r}{n_1,n_2,\ldots,n_r}=\frac{(n_1+n_2+\cdots+n_r)!}{n_1!\,n_2!\,\cdots\, n_r!}.$$

In the question this simplifies to

$$P(4,4,\ldots,4) = \frac{52!}{(4!)^{13}}.$$

The desired probability is obtained by dividing the number of derangements by the number of (equally likely) permutations,

$$\frac{D(4,4,\ldots,4)}{P(4,4,\ldots,4)} = \left|\frac{52! + p_{51}51! + \cdots + p_1 1! + p_0 0!}{52!}\right|.$$

I computed this with Wolfram Alpha.

Answer 2

You can do this with a binomial distribution. I used R.

Every time you turn over a card, you have a 1 in 13 chance of it matching the number you count out. You need to turn over 52 cards in a row, without a match. (You can think of it as four sets of 13, that doesn't change it).

In R (or whatever you use probably has a similar function).

> pbinom(q = 0, size = 52, prob = 1/13)
[1] 0.01557294

Equivalently, you need an event that has a 12/13 chance of happening (no match), but you need it to happen 52 times in a row.

> (12/13)^52
[1] 0.01557294

Very close, but not quite the same as you got in your simulation.

(Hmmm... my simulation gives a slightly higher number, not as high as yours. I'm not sure why.

> library(dplyr)
> CountOutCards <- function (x) {
+   sum(rep(1:13, 4) == sample(rep(1:13, 4), 52, FALSE))
+ }
> 
> x <- (lapply(1:1e6, CountOutCards) ) %>% unlist 
> mean(x == 0)
[1] 0.016093

)