What are the odds of winning this solitaire card game?

Question

I remember a certain very simple solitaire game played on a standard 52 card deck that goes as follows: Shuffle the deck, place the top card face up, and call out "ace." If the card is an ace, you lose. If not, proceed with the next card and call out "two." If the card is a two, you lose. If not, proceed calling out 3, 4, 5, 6, 7, 8, 9, 10, J, Q, K, then A, 2, 3 and so on. You'll cycle through the ranks in order exactly four times (if you don't lose), and suit doesn't matter. The object is to exhaust the deck without ever calling out the rank of the card you're exposing. The question is, what are the odds of winning this game?

First card is easy. There's a $\frac{4}{52}=\frac{1}{13}$ chance of exposing an ace, so a $\frac{12}{13}$ chance of getting past the first card. But now there's either a $\frac{3}{51}$ or $\frac{4}{51}$ chance of exposing a two, depending on whether the first card exposed was a two, so: $$\frac{1}{13}\left(\frac{3}{51}\right)+\frac{12}{13}\left(\frac{4}{51}\right)$$ Well this just gives $\frac{1}{13}$ again, so a $\frac{12}{13}$ chance of moving on. Could it be that at any draw, the odds are $\frac{1}{13}$ of calling out a match? Intuitively I'm struggling with this, but going on this assumption, I'm coming up with the odds of winning this game to be $$\left(\frac{12}{13}\right)^{52}\approx0.0156$$ or roughly 1.5%.

Those odds are more or less consistent with my experience playing this game, but are they correct?

EDIT: A solution from another post has been offered which appears correct, but I'd like to expand my question to ask the odds of drawing/not drawing a match at any given iteration (assuming you haven't yet lost). For example, are my odds for the second iteration correct? And what does the sequence of odds look like?

EDIT 2: Calculating the odds that the game ends with the n-th draw (or $D_n$), I came up with the following for the first three draws. Draw 1 is obvious: $$P(D_1)=\frac{4}{52}=\frac{1}{13}$$ Draw 2 I've already established above except I corrected the error that @saulspatz pointed out, and for clarity I've stopped simplifying: $$P(D_2)=\frac{4}{52}\left(\frac{3}{51}\right)+\frac{44}{52}\left(\frac{4}{51}\right)$$ For Draw 3, there are five ways the first two draws can play out (without losing) with odds shown below, and 50 cards remain after each scenario: $$D_1=3;\quad D_2=3 \qquad \frac{4}{52}\cdot\frac{3}{51}$$ $$D_1=3;\quad D_2\neq2,3 \qquad \frac{4}{52}\cdot\frac{44}{51}$$ $$D_1\neq A,3;\quad D_2=3 \qquad \frac{44}{52}\cdot\frac{4}{51}$$ $$D_1=2;\quad D_2\neq 2,3 \qquad \frac{4}{52}\cdot\frac{44}{51}$$ $$D_1\neq A,2,3;\quad D_2\neq 2,3 \qquad \frac{40}{52}\cdot\frac{43}{51}$$

So after Scenario 1, two 3's remain; after Scenario's 2 and 3, three 3's remain; and after Scenario's 4 and 5, four 3's remain. And so: $$P(D_3)= \frac{4}{52}\cdot\frac{3}{51} \left(\frac{2}{50}\right)+ \frac{4}{52}\cdot\frac{44}{51} \left(\frac{3}{50}\right)+ \frac{44}{52}\cdot\frac{4}{51} \left(\frac{3}{50}\right)+ \frac{4}{52}\cdot\frac{44}{51} \left(\frac{4}{50}\right)+ \frac{40}{52}\cdot\frac{43}{51} \left(\frac{4}{50}\right) $$ Not sure if or how much this helps, but I felt the need to lay it all out. Anything beyond $D_3$ gets unwieldy pretty fast. I still feel like they're must be some reasonably compact formula that shows the odds of losing at $D_n$ but I'm not equipped to derive it. But it looks like there might be some ways to generalize what's happening from looking at the first three formulas.

Answer 1

Your calculation of the probability of getting a $2$ on the second card is erroneous. You are right in saying that with probability $\frac1{13}$, the first card was a $2$, and in that case the probability of drawing a $2$ is $\frac3{51}$ and also that if you didn't draw a $2$ on the first card, the probability of drawing one on the second is $\frac4{51}$. However, you should not be weighting this by $\frac{12}{13}$ but by $\frac{11}{13}$, since with probability $\frac1{13}$ you draw an Ace on the first card and never get to the second card.

The general problem you raise, which we may state as finding the probability that the game ends on card $n$ for $1\leq n\leq52$ seems difficult by straightforward calculation. If we want to compute the probability that the $17$ card is a $4$, we must know the probabilities that there are $4,3,2,1$, or $0$ fours left in the deck, and each of these depends on how many of each other rank there are. Basically, for every choice of the first $16$ cards, we have to compute the probability of getting past them all.

The problem can be solved with rook polynomials, though. For the example of the $17$th card, we are considering first hitting on the second $4$. For the first $16$ spots, we have the usual prohibitions. For the $17$th spot, we prohibit any card but a $4$. For the $18$th and later spots, there are no prohibitions.

It turns out to be rather easy to write down the formulas for the appropriate rook polynomials. There are two main facts used to compute rook polynomials, as explained here. If the board splits into disjunct sub-boards, then the rook polynomial of the board equals the product of the rook polynomials of the sub-boards. (Two sub-boards are disjunct if if no cell of one of them is in the same row or column as a cell in the other.)

The second fact, explained in the reference by the equation $$R(x,B)=R(x,B_s)+xR(x,B_s^*)\tag1$$ is that if $B$ is a board and $s$ a cell of $B$, then the rook polynomial of $B$ is the rook polynomial of $B$ with $s$ deleted ($B_s$) plus $x$ times the rook polynomial of $B$ with both the row and columns of $s$ deleted ($B_s^*$). The terms correspond to the cases where a rook is not placed on $s$ and where one is so placed.

In this problem, the boards will consist largely of disjunct rectangles so the other fact we need is that the rook polynomial of an $m\times n$ rectangle is $$\sum_{k=0}^{\min{m,n}}\binom mk\binom nk k!x^k$$ Clearly, we can place at most $\min(n,m)$ non-attacking rooks on an $m\times n$ rectangle. To place $k$ rooks, there are $\binom mk$ ways to choose their rows, and $\binom nk$ ways to choose their columns. The selected rows and columns intersect in a $k\times k$ grid, and there are $k!$ ways to place $k$ non-attacking rooks on such a grid.

Now let's look at the example of the $17$th card again. After permuting rows and columns, we have three $2\times4$ rectangles, corresponding to Aces, deuces and treys, and ten $1\times4$ rectangles, for the other ranks. These are pairwise disjunct. We have one row corresponding to the second four, and it has $48$ cells. The rectangle corresponding to the first four is disjunction form this row also, so by the first rule, we can pull out the columns corresponding to the first four, and multiply by $1+4x$, the rook polynomial of a $1\times4$ rectangle.

Now we repeatedly apply formula $(1)$ removing the squares of the $1\times 48$ rectangle one after another. When we do this, the second term in the formula is easy to compute, because when we delete the last row, the board is a union of disjunct rectangles, and we can compute the rook polynomial immediately.

Twelve times we are left with nine $1\times4$ rectangles, two $2\times4$ rectangles, and one $2\times3$ rectangles. The rook polynomial of such a board is $$p_1(x):=\big(1+4x\big)^9\left(1+8x+12x^2\right)^2\left(1+6x+6x^2\right)$$

The other $36$ times, we are left with three $2\times4$ rectangles, eight $1\times4$ rectangles, and one $1\times3$ rectangle. The rook polynomial of such a board is $$p_2(x):=\big(1+4x\big)^8\left(1+8x+12x^2\right)^3\big(1+3x\big)$$

When we have deleted each of the $48$ cells in the row, we are left with a board with three $2\times4$ rectangles and $9$ $1\times4$ rectangles and a rook polynomial of $$p_3(x):=\big(1+4x\big)^9\left(1+8x+12x^2\right)^3$$

Putting it all together, the rook polynomial of the board is $$(1+4x)(12xp_1(x)+36xp_2(x)+p_3(x))$$

To actually compute the coefficients and compute the probabilities, a computer program is desirable.

I'm still planning to do this, but I'm not sure I'll get to it today, so I'm posting this in the meantime.

EDIT

I've got my script working, though it took an embarrassingly long time to find the bug. I've tried to comment it very thoroughly. The script counts the number of occurrences of each case (we successfully get past exactly $n$ cards for $0\leq n\leq 52$ and checks that they add up to $52!$, so I'm confident that it's correct.

Here's the script:

'''
Rencontre played with a deck of s suits and r ranks.
What is the probability that exactly n cards are unmatched
before the game ends, 0 <= n <= rs?
'''
from sympy import binomial, Poly, symbols, factorial

x = symbols('x')

def rectangle(n,m):
    '''
    Rook polynomial of n x m rectangle
    '''
    if min(n,m)==0:
        return 1
    return sum([binomial(n,k)*binomial(m,k)*factorial(k)*x**k for k in range(1+min(n,m)) ]).as_poly()

class RookPoly:
    def __init__(self,ranks=13, suits=4):
        self. ranks = ranks
        self.suits = suits
        self.cards = ranks * suits
        self.deals = factorial(self.cards)
        
    def rp(self, n):
        r = self.ranks
        s = self.suits        
        if n == 0:
            return rectangle(1,s*(r-1))
        
        if n == self.cards:
            return rectangle(s,s)**r
    
        a, b = divmod(n, r)
        # We get past b ranks a+1 times, and
        # we get past r-b ranks a times.
        
        if a == 0:
            # We don't get past all the ranks
            # We have b 1 x s disjunct rectangles
            # and the long row has (r-1)*s cells.
            # (r-b-1)s times,  when we delete a cell from the
            # long row, we delete no cells in any other row.
            # We also get this board when we are done with
            # the long row
        
            p1 = (x*(r-b-1)*s+1)*rectangle(1,s)**b
    
            # sb times we get 
            #    - b-1 1 x s rectangles
            #     - one 1x s-1 rectangle
            
            p2 = b*s*x*rectangle(1,s)**(b-1)*rectangle(1,s-1) 
            return p1+p2
        
        # Each rank has been successfully passed.  The rectangle corresponsing
        # to the first card hit is disjunct from the rest of the board.
        
        p = rectangle(a,s)
        
        # The long row has s*(r-1) cells.
        # When we delete one of it cells
        # if b > 0, then b*s times we leave 
        #     - b-1 (a+1) x s rectangles
        #     - one (a+1) x (s-1) rectangle
        #     - r-b-1 a x s rectangles 
        
        p1 =  b*s*rectangle(a+1,s)**(b-1)*rectangle(a+1,s-1)*rectangle(a,s)**(r-b-1) if b > 0 else 0
        
        # if b < r-1, then
        # (r-1-b)*s times we leave
        #      - b (a+1) x s rectangles
        #      - r-b-2 a x s rectangles
        #      - one a x (s-1) rectangle
        
        p2 = (r-1-b)*s*rectangle(a+1,s)**b*rectangle(a,s)**(r-b-2)*rectangle(a,s-1) if b < r-1 else 0
        
        # When the long row is eliminated, we have
        #       - b (a+1) x s rectangles 
        #        - r-b-1 a x s rectangles
        
        p3 = rectangle(a+1,s)**b*rectangle(a,s)**(r-b-1)
        
        return p*(p1*x+p2*x+p3)
    
    def occurs(self, n):
        '''
        Number of deals in which we successfully pass
        n cards, 0 <= n <= self.cards
        '''
        rook = self.rp(n)
        
        # Inclusion and exclusion
        n = self.cards
        coeffs = rook.all_coeffs()[::-1]
        return sum((-1)**k*coeffs[k]*factorial(n-k) for k in range(1+rook.degree()))
    
    def prob(self, n):
        successes = self.occurs(n)
        return successes, successes/self.deals
    
    def allProbs(self):
        results = [self.prob(n) for n in range(1+self.cards)]
        cases = sum(result[0] for result in results)
        try:
            assert cases == self.deals
            print(f'{self.ranks} ranks, {self.suits} suits')
            print('Passed audit')
        except AssertionError:
            print('FAILED AUDIT')
            print(f'{cases} cases, {self.deals} deals')
        return results

R = RookPoly(ranks=13, suits=4)
results = R.allProbs()
for n in range(53):
    print('%2d %.9f'%(n, results[n][1]))
    

and here's the output:

13 ranks, 4 suits
Passed audit
 0 0.076923077
 1 0.070889894
 2 0.065339367
 3 0.060232093
 4 0.055531956
 5 0.051205844
 6 0.047223392
 7 0.043556750
 8 0.040180375
 9 0.037070829
10 0.034206605
11 0.031567966
12 0.029136795
13 0.027456543
14 0.025303032
15 0.023321812
16 0.021498817
17 0.019821152
18 0.018276998
19 0.016855515
20 0.015546761
21 0.014341617
22 0.013231716
23 0.012209382
24 0.011267570
25 0.010399811
26 0.009800333
27 0.009031639
28 0.008324448
29 0.007673738
30 0.007074905
31 0.006523730
32 0.006016344
33 0.005549197
34 0.005119034
35 0.004722868
36 0.004357960
37 0.004021793
38 0.003712059
39 0.003498176
40 0.003223787
41 0.002971354
42 0.002739082
43 0.002525330
44 0.002328589
45 0.002147480
46 0.001980734
47 0.001827190
48 0.001685781
49 0.001555530
50 0.001435539
51 0.001324982
52 0.016232727

It's interesting that $\Pr(N=n)$ is a strictly decreasing function of $n$ until it takes a sudden leap upward at $n=52$. It is know that if the deck has $r$ ranks, and $s$ suits, and let $r$ increase while holding $s$ fixed, then the probability of getting through the whole deck goes to $e^{-s}$ as $r\to\infty$. I wonder what, if anything, can be said about the other probabilities.

I really enjoyed this. Thank you for the problem.