Expected number I will be on after drawing cards until I get an ace, 2, 3, and so forth

Question

I am having some trouble solving the following.

You draw cards from a standard 52-card deck without replacement until you get an ace. You draw from what is remaining until you get a 2. You continue on with 3. What is the expected number you will be on after the entire deck runs out?

It was natural to let

• $T_i = \text{first position of card whose value is }i$
• $U_i = \text{last position of card whose value is }i$

So the problem essentially amounts to figuring out the probability that you will be on $k$ when the deck runs out, namely:

$$Pr(T_1<\cdots<T_k \cap U_{k+1} < T_k)$$

I can see that

$$Pr(T_1<\cdots<T_k) = 1/k! \\ \text{and} \\ Pr(U_{k+1} < T_k) = 1/70$$

but could not get any further...

Answer 1

following @gung's idea, i believe the expected value would be 5.84? and from my interpretation of the comments, i'm assuming "A" is an almost impossible value (unless the last four cards in the deck are all aces). here are the results of a 100,000 iteration monte carlo simulation

results
    2     3     4     5     6     7     8     9     J     K     Q     T 
 1406  7740 16309 21241 19998 15127  9393  4906   976   190   380  2334 

and here's the R code in case you'd like to play with it..

# monte carlo card-drawing functions from here
# http://streaming.stat.iastate.edu/workshops/r-intro/lectures/5-Rprogramming.pdf

# create a straightforward deck of cards
create_deck <-
    function( ){
        suit <- c( "H" , "C" , "D" , "S" )
        rank <- c( "A" , 2:9 , "T" , "J" , "Q" , "K" )
        deck <- NULL
        for ( r in rank ) deck <- c( deck , paste( r , suit ) )
        deck
    }

# construct a function to shuffle everything
shuffle <- function( deck ){ sample( deck , length( deck ) ) }

# draw one card at a time
draw_cards <-
    function( deck , start , n = 1 ){
        cards <- NULL

        for ( i in start:( start + n - 1 ) ){
            if ( i <= length( deck ) ){
                cards <- c( cards , deck[ i ] )
            }
        }

        return( cards )
    }

# create an empty vector for your results
results <- NULL

# run your simulation this many times..
for ( i in seq( 100000 ) ){
    # create a new deck
    sdeck <- shuffle( create_deck() )

    d <- sdeck[ grep('A|2' , sdeck ) ]
    e <- identical( grep( "2" , d ) , 1:4 )

    # loop through ranks in this order
    rank <- c( "A" , 2:9 , "T" , "J" , "Q" , "K" )

    # start at this position
    card.position <- 0

    # start with a blank current.draw
    current.draw <- ""

    # start with a blank current rank
    this.rank <- NULL

    # start with the first rank
    rank.position <- 1

    # keep drawing until you find the rank you wanted
    while( card.position < 52 ){

        # increase the position by one every time
        card.position <- card.position + 1

        # store the current draw for testing next time
        current.draw <- draw_cards( sdeck , card.position )

        # if you draw the current rank, move to the next.
        if ( grepl( rank[ rank.position ] , current.draw ) ) rank.position <- rank.position + 1

        # if you have gone through every rank and are still not out of cards,
        # should it still be a king?  this assumes yes.
        if ( rank.position == length( rank ) ) break        

    }

    # store the rank for this iteration.
    this.rank <- rank[ rank.position ]

    # at the end of the iteration, store the result
    results <- c( results , this.rank )

}

# print the final results
table( results )

# make A, T, J, Q, K numerics
results[ results == 'A' ] <- 1
results[ results == 'T' ] <- 10
results[ results == 'J' ] <- 11
results[ results == 'Q' ] <- 12
results[ results == 'K' ] <- 13
results <- as.numeric( results )

# and here's your expected value after 100,000 simulations.
mean( results )

Answer 2

For a simulation it's crucial to be correct as well as fast. Both these objectives suggest writing code that targets core capabilities of the programming environment as well as code that is as short and simple as possible, because simplicity lends clarity and clarity promotes correctness. Here is my attempt to achieve both in R:

#
# Simulate one play with a deck of `n` distinct cards in `k` suits.
#
sim <- function(n=13, k=4) {
  deck <- sample(rep(1:n, k)) # Shuffle the deck
  deck <- c(deck, 1:n)        # Add sentinels to terminate the loop
  k <- 0                      # Count the cards searched for
  for (j in 1:n) {
    k <- k+1                          # Count this card
    deck <- deck[-(1:match(j, deck))] # Deal cards until `j` is found
    if (length(deck) < n) break       # Stop when sentinels are reached
  }
  return(k)                   # Return the number of cards searched
}

Applying this in a reproducible way can be done with the replicate function after setting the random number seed, as in

> set.seed(17);  system.time(d <- replicate(10^5, sim(13, 4)))
   user  system elapsed 
   5.46    0.00    5.46

That's slow, but fast enough to conduct fairly lengthy (and therefore precise) simulations repeatedly without waiting. There are several ways we can exhibit the result. Let's start with its mean:

> n <- length(d)
> mean(d)
[1] 5.83488

> sd(d) / sqrt(n)
[1] 0.005978956

The latter is the standard error: we expect the simulated mean to be within two or three SEs of the true value. That places the true expectation somewhere between $5.817$ and $5.853$.

We might also want to see a tabulation of the frequencies (and their standard errors). The following code prettifies the tabulation a little:

u <- table(d)
u.se <- sqrt(u/n * (1-u/n)) / sqrt(n)
cards <- c("A", "2", "3", "4", "5", "6", "7", "8", "9", "T", "J", "Q", "K")
dimnames(u) <- list(sapply(dimnames(u), function(x) cards[as.integer(x)]))
print(rbind(frequency=u/n, SE=u.se), digits=2)

Here is the output:

                2       3      4      5      6      7       8       9       T       J       Q       K
frequency 0.01453 0.07795 0.1637 0.2104 0.1995 0.1509 0.09534 0.04995 0.02249 0.01009 0.00345 0.00173
SE        0.00038 0.00085 0.0012 0.0013 0.0013 0.0011 0.00093 0.00069 0.00047 0.00032 0.00019 0.00013

---

How can we know the simulation is even correct? One way is to test it exhaustively for smaller problems. For that reason this code was written to attack a small generalization of the problem, replacing $13$ distinct cards with n and $4$ suits with k. However, for the testing it is important to be able to feed the code a deck in a predetermined order. Let's write a slightly different interface to the same algorithm:

draw <- function(deck) {
  n <- length(sentinels <- sort(unique(deck)))
  deck <- c(deck, sentinels)
  k <- 0
  for (j in sentinels) {
    k <- k+1
    deck <- deck[-(1:match(j, deck))]
    if (length(deck) < n) break
  }
  return(k)
}

(It is possible to use draw in place of sim everywhere, but the extra work done at the beginning of draw makes it twice as slow as sim.)

We can use this by applying it to every distinct shuffle of a given deck. Since the purpose here is just a few one-off tests, efficiency in generating those shuffles is unimportant. Here is a quick brute-force way:

n <- 4 # Distinct cards
k <- 2 # Number of suits
d <- expand.grid(lapply(1:(n*k), function(i) 1:n))
e <- apply(d, 1, function(x) var(tabulate(x))==0)
g <- apply(d, 1, function(x) length(unique(x))==n)
d <- d[e & g,]

Now d is a data frame whose rows contain all the shuffles. Apply draw to each row and count the results:

d$result <- apply(as.matrix(d), 1, draw)
    (counts <- table(d$result))

The output (which we will use in a formal test momentarily) is

   2    3    4 
 420  784 1316 

(The value of $420$ is easy to understand, by the way: we would still be working on card $2$ if and only if all the twos preceded all the aces. The chance of this happening (with two suits) is $1/\binom{2+2}{2} = 1/6$. Out of the $2520$ distinct shuffles, $2520/6 = 420$ have this property.)

We can test the output with a chi-squared test. To this end I apply sim $10,000$ times to this case of $n = 4$ distinct cards in $k = 2$ suits:

>set.seed(17)
>d.sim <- replicate(10^4, sim(n, k))
>print((rbind(table(d.sim) / length(d.sim), counts / dim(d)[1])), digits=3)

         2     3     4
[1,] 0.168 0.312 0.520
[2,] 0.167 0.311 0.522

> chisq.test(table(d.sim), p=counts / dim(d)[1])

    Chi-squared test for given probabilities

data:  table(d.sim) 
X-squared = 0.2129, df = 2, p-value = 0.899

Because $p$ is so high, we find no significant difference between what sim says and the values computed by exhaustive enumeration. Repeating this exercise for some other (small) values of $n$ and $k$ produces comparable results, giving us ample reason to trust sim when applied to $n=13$ and $k=4$.

Finally, a two-sample chi-squared test will compare the output of sim to the output reported in another answer:

>y <- c(1660,8414,16973,21495,20021,14549,8957,4546,2087,828,313,109)
>chisq.test(cbind(u, y))

data:  cbind(u, y) 
X-squared = 142.2489, df = 11, p-value < 2.2e-16

The enormous chi-squared statistic produces a p-value that is essentially zero: without a doubt, sim disagrees with the other answer. There are two possible resolutions of the disagreement: one (or both!) of these answers is incorrect or they implement different interpretations of the question. For instance, I have interpreted "after the deck runs out" to mean after observing the last card and, if allowable, updating the "number you will be on" before terminating the procedure. Conceivably that last step was not meant to be taken. Perhaps some such subtle difference of interpretation will explain the disagreement, at which point we can modify the question to make it clearer what is being asked.

Answer 3

There is an exact answer (in the form of a matrix product, presented in point 4 below). A reasonably efficient algorithm to compute it exists, deriving from these observations:

1. A random shuffle of $N+k$ cards can be generated by randomly shuffling $N$ cards and then randomly interspersing the remaining $k$ cards within them.

2. By shuffling only the aces, and then (applying the first observation) interspersing the twos, then the threes, and so on, this problem can be viewed as a chain of thirteen steps.

3. We do need to keep track of more than the value of the card we are seeking. When doing this, however, we don't need to account for the position of the mark relative to all the cards, but only its position relative to cards of equal or smaller value.

   Imagine placing a mark on the first ace, and then marking the first two found after it, and so on. (If at any stage the deck runs out without displaying the card we are currently seeking, we will leave all cards unmarked.) Let the "place" of each mark (when it exists) be the number of cards of equal or lower value that were dealt when the mark was made (including the marked card itself). The places contain all the essential information.

4. The place after the $i^\text{th}$ mark is made is a random number. For a given deck, the sequence of these places forms a stochastic process. It in fact is a Markov process (with variable transition matrix). An exact answer can therefore be calculated from twelve matrix multiplications.

Using these ideas, this machine obtains a value of $5.8325885529019965$ (computing in double precision floating point) in $1/9$ second. This approximation of the exact value $$\frac{1982600579265894785026945331968939023522542569}{339917784579447928182134345929899510000000000}$$ is accurate to all digits shown.

The rest of this post provides details, presents a working implementation (in R), and concludes with some comments about the question and the efficiency of the solution.

---

Generating random shuffles of a deck

It is actually clearer conceptually and no more complicated mathematically to consider a "deck" (aka multiset) of $N = k_1+k_2+\cdots+k_m$ cards of which there are $k_1$ of the lowest denomination, $k_2$ of the next lowest, and so on. (The question as asked concerns the deck determined by the $13$-vector $(4,4,\ldots,4)$.)

A "random shuffle" of $N$ cards is one permutation taken uniformly and randomly from the $N! = N\times(N-1)\times\cdots\times 2\times 1$ permutations of the $N$ cards. These shuffles fall into groups of equivalent configurations because permuting the $k_1$ "aces" among themselves changes nothing, permuting the $k_2$ "twos" among themselves also changes nothing, and so on. Therefore each group of permutations that look identical when the suits of the cards are ignored contains $k_1!\times k_2!\times \cdots \times k_m!$ permutations. These groups, whose number therefore is given by the multinomial coefficient

$$\binom{N}{k_1,k_2,\ldots,k_m} = \frac{N!}{k_1!k_2!\cdots k_m!},$$

are called "combinations" of the deck.

There is another way to count the combinations. The first $k_1$ cards can form only $k_1!/k_1! = 1$ combination. They leave $k_1+1$ "slots" between and around them in which the next $k_2$ cards can be placed. We could indicate this with a diagram where "$*$" designates one of the $k_1$ cards and "$\_$" designates a slot that can hold between $0$ and $k_2$ additional cards:

$$\underbrace{\_*\_*\_\cdots\_*\_}_{k_1\text{ stars}}$$

When $k_2$ additional cards are interspersed, the pattern of stars and new cards partitions the $k_1+k_2$ cards into two subsets. The number of distinct such subsets is $\binom{k_1+k_2}{k_1,k_2} = \frac{(k_1+k_2)!}{k_1!k_2!}$.

Repeating this procedure with $k_3$ "threes," we find there are $\binom{(k_1+k_2)+k_3}{k_1+k_2,k_3}= \frac{(k_1+k_2+k_3)!}{(k_1+k_2)!k_3!}$ ways to intersperse them among the first $k_1+k_2$ cards. Therefore the total number of distinct ways to arrange the first $k_1+k_2+k_3$ cards in this manner equals

$$1\times\frac{(k_1+k_2)!}{k_1!k_2!}\times\frac{(k_1+k_2+k_3)!}{(k_1+k_2)!k_3!} = \frac{(k_1+k_2+k_3)!}{k_1!k_2!k_3!}.$$

After finishing the last $k_n$ cards and continuing to multiply these telescoping fractions, we find that the number of distinct combinations obtained equals the total number of combinations as previously counted, $\binom{N}{k_1,k_2,\ldots,k_m}$. Therefore we have overlooked no combinations. That means this sequential process of shuffling the cards correctly captures the probabilities of each combination, assuming that at each stage each possible distinct way of interspersing the new cards among the old is taken with uniformly equal probability.

The place process

Initially, there are $k_1$ aces and obviously the very first is marked. At later stages there are $n = k_1 + k_2 + \cdots + k_{j-1}$ cards, the place (if a marked card exists) equals $p$ (some value from $1$ through $n$), and we are about to intersperse $k=k_j$ cards around them. We can visualize this with a diagram like

$$\underbrace{\_*\_*\_\cdots\_*\_}_{p-1\text{ stars}}\odot\underbrace{\_*\_\cdots\_*\_}_{n-p\text{ stars}}$$

where "$\odot$" designates the currently marked symbol. Conditional on this value of the place $p$, we wish to find the probability that the next place will equal $q$ (some value from $1$ through $n+k$; by the rules of the game, the next place must come after $p$, whence $q\ge p+1$). If we can find how many ways there are to intersperse the $k$ new cards in the blanks so that the next place equals $q$, then we can divide by the total number of ways to intersperse these cards (equal to $\binom{n+k}{k}$, as we have seen) to obtain the transition probability that the place changes from $p$ to $q$. (There will also be a transition probability for the place to disappear altogether when none of the new cards follow the marked card, but there is no need to compute this explicitly.)

Let's update the diagram to reflect this situation:

$$\underbrace{\_*\_*\_\cdots\_*\_}_{p-1\text{ stars}}\odot\underbrace{**\cdots*}_{s\text{ stars}}\ \vert\ \underbrace{\_*\_\cdots\_*\_}_{n-p-s\text{ stars}}$$

The vertical bar "$\vert$" shows where the first new card occurs after the marked card: no new cards may therefore appear between the $\odot$ and the $\vert$ (and therefore no slots are shown in that interval). We do not know how many stars there are in this interval, so I have just called it $s$ (which may be zero) The unknown $s$ will disappear once we find the relationship between it and $q$.

Suppose, then, we intersperse $j$ new cards around the stars before the $\odot$ and then--independently of that--we intersperse the remaining $k-j-1$ new cards around the stars after the $\vert$. There are

$$\tau_{n,k}(s,p) = \binom{(p-1)+j}{j}\binom{(n-p-s) + (k-j)-1}{k-j-1}$$

ways to do this. Notice, though--this is the trickiest part of the analysis--that the place of $\vert$ equals $p+s+j+1$ because

• There are $p$ "old" cards at or before the mark.
• There are $s$ old cards after the mark but before $\vert$.
• There are $j$ new cards before the mark.
• There is the new card represented by $\vert$ itself.

Thus, $\tau_{n,k}(s,p)$ gives us information about the transition from place $p$ to place $q=p+s+j+1$. When we track this information carefully for all possible values of $s$, and sum over all these (disjoint) possibilities, we obtain the conditional probability of place $q$ following place $p$,

$${\Pr}_{n,k}(q|p) = \left(\sum_j \binom{p-1+j}{j}\binom{n+k-q}{k-j-1}\right) / \binom{n+k}{k}$$

where the sum starts at $j=\max(0, q-(n+1))$ and ends at $j=\min(k-1, q-(p+1)$. (The variable length of this sum suggests there is unlikely to be a closed formula for it as a function of $n, k, q,$ and $p$, except in special cases.)

The algorithm

Initially there is probability $1$ that the place will be $1$ and probability $0$ it will have any other possible value in $2, 3, \ldots, k_1$. This can be represented by a vector $p_1 = (1, 0, \ldots, 0)$.

After interspersing the next $k_2$ cards, the vector $p_1$ is updated to $p_2$ by multiplying it (on the left) by the transition matrix $(\Pr_{k_1,k_2}(q|p), 1\le p\le k_1, 1\le q\le k_2)$. This is repeated until all $k_1+k_2+\cdots+k_m$ cards have been placed. At each stage $j$, the sum of the entries in the probability vector $p_j$ is the chance that some card has been marked. Whatever remains to make the value equal to $1$ therefore is the chance that no card is left marked after step $j$. The successive differences in these values therefore give us the probability that we could not find a card of type $j$ to mark: that is the probability distribution of the value of the card we were looking for when the deck runs out at the end of the game.

---

Implementation

The following R code implements the algorithm. It parallels the preceding discussion. First, calculation of the transition probabilities is performed by t.matrix (without normalization with the division by $\binom{n+k}{k}$, making it easier to track the calculations when testing the code):

t.matrix <- function(q, p, n, k) {
  j <- max(0, q-(n+1)):min(k-1, q-(p+1))
  return (sum(choose(p-1+j,j) * choose(n+k-q, k-1-j))
}

This is used by transition to update $p_{j-1}$ to $p_j$. It calculates the transition matrix and performs the multiplication. It also takes care of computing the initial vector $p_1$ if the argument p is an empty vector:

#
# `p` is the place distribution: p[i] is the chance the place is `i`.
#
transition <- function(p, k) {
  n <- length(p)
  if (n==0) {
    q <- c(1, rep(0, k-1))
  } else {
    #
    # Construct the transition matrix.
    #
    t.mat <- matrix(0, nrow=n, ncol=(n+k))
    #dimnames(t.mat) <- list(p=1:n, q=1:(n+k))
    for (i in 1:n) {
      t.mat[i, ] <- c(rep(0, i), sapply((i+1):(n+k), 
                                        function(q) t.matrix(q, i, n, k)))
    }
    #
    # Normalize and apply the transition matrix.
    #
    q <- as.vector(p %*% t.mat / choose(n+k, k))
  }
  names(q) <- 1:(n+k)
  return (q)
}

We can now easily compute the non-mark probabilities at each stage for any deck:

#
# `k` is an array giving the numbers of each card in order;
# e.g., k = rep(4, 13) for a standard deck.
#
# NB: the *complements* of the p-vectors are output.
#
game <- function(k) {
  p <- numeric(0)
  q <- sapply(k, function(i) 1 - sum(p <<- transition(p, i)))
  names(q) <- names(k)
  return (q)
}

Here they are for the standard deck:

k <- rep(4, 13)
names(k) <- c("A", 2:9, "T", "J", "Q", "K")
(g <- game(k))

The output is

         A          2          3          4          5          6          7          8          9          T          J          Q          K 
0.00000000 0.01428571 0.09232323 0.25595013 0.46786622 0.66819134 0.81821790 0.91160622 0.96146102 0.98479430 0.99452614 0.99818922 0.99944610

According to the rules, if a king was marked then we would not look for any further cards: this means the value of $0.9994461$ has to be increased to $1$. Upon doing so, the differences give the distribution of the "number you will be on when the deck runs out":

> g[13] <- 1; diff(g)
          2           3           4           5           6           7           8           9           T           J           Q           K 
0.014285714 0.078037518 0.163626897 0.211916093 0.200325120 0.150026562 0.093388313 0.049854807 0.023333275 0.009731843 0.003663077 0.001810781

(Compare this to the output I report in a separate answer describing a Monte-Carlo simulation: they appear to be the same, up to expected amounts of random variation.)

The expected value is immediate:

> sum(diff(g) * 2:13)
[1] 5.832589

All told, this required only a dozen lines or so of executable code. I have checked it against hand calculations for small values of $k$ (up to $3$). Thus, if any discrepancy becomes apparent between the code and the preceding analysis of the problem, trust the code (because the analysis may have typographical errors).

---

Remarks

Relationships to other sequences

When there is one of each card, the distribution is a sequence of reciprocals of whole numbers:

> 1/diff(game(rep(1,10)))
[1]      2      3      8     30    144    840   5760  45360 403200

The value at place $i$ is $i! + (i-1)!$ (starting at place $i=1$). This is sequence A001048 in the Online Encyclopedia of Integer Sequences. Accordingly, we might hope for a closed formula for the decks with constant $k_i$ (the "suited" decks) that would generalize this sequence, which itself has some profound meanings. (For instance, it counts sizes of the largest conjugacy classes in permutation groups and is also related to trinomial coefficients.) (Unfortunately, the reciprocals in the generalization for $k\gt 1$ are not usually integers.)

The game as a stochastic process

Our analysis makes it clear that the initial $i$ coefficients of the vectors $p_j$, $j\ge i$, are constant. For example, let's track the output of game as it processes each group of cards:

> sapply(1:13, function(i) game(rep(4,i)))

[[1]]
[1] 0

[[2]]
[1] 0.00000000 0.01428571

[[3]]
[1] 0.00000000 0.01428571 0.09232323

[[4]]
[1] 0.00000000 0.01428571 0.09232323 0.25595013

...

[[13]]
 [1] 0.00000000 0.01428571 0.09232323 0.25595013 0.46786622 0.66819134 0.81821790 0.91160622 0.96146102 0.98479430 0.99452614 0.99818922 0.99944610

For instance, the second value of the final vector (describing the results with a full deck of 52 cards) already appeared after the second group was processed (and equals $1/\binom{8}{4}=1/70$). Thus, if you want information only about the marks up through the $j^\text{th}$ card value, you only have to perform the calculation for a deck of $k_1+k_2+\cdots+k_j$ cards.

Because the chance of not marking a card of value $j$ is getting quickly close to $1$ as $j$ increases, after $13$ types of cards in four suits we have almost reached a limiting value for the expectation. Indeed, the limiting value is approximately $5.833355$ (computed for a deck of $4 \times 32$ cards, at which point double precision rounding error prevents going any further).

Timing

Looking at the algorithm applied to the $m$-vector $(k,k, \ldots, k)$, we see its timing should be proportional to $k^2$ and--using a crude upper bound--not any worse than proportional to $m^3$. By timing all calculations for $k=1$ through $7$ and $n=10$ through $30$, and analyzing only those taking relatively long times ($1/2$ second or longer), I estimate the computation time is approximately $O(k^2 n^{2.9})$, supporting this upper-bound assessment.

One use of these asymptotics is to project calculation times for larger problems. For instance, seeing that the case $k=4, n=30$ takes about $1.31$ seconds, we would estimate that the (very interesting) case $k=1, n=100$ would take about $1.31(1/4)^2(100/30)^{2.9}\approx 2.7$ seconds. (It actually takes $2.87$ seconds.)

Answer 4

Hacked a simple Monte Carlo in Perl and found approximately $5.8329$.

#!/usr/bin/perl

use strict;

my @deck = (1..13) x 4;

my $N = 100000; # Monte Carlo iterations.

my $mean = 0;

for (my $i = 1; $i <= $N; $i++) {
    my @d = @deck;
    fisher_yates_shuffle(\@d);
    my $last = 0;
        foreach my $c (@d) {
        if ($c == $last + 1) { $last = $c }
    }
    $mean += ($last + 1) / $N;
}

print $mean, "\n";

sub fisher_yates_shuffle {
    my $array = shift;
        my $i = @$array;
        while (--$i) {
        my $j = int rand($i + 1);
        @$array[$i, $j] = @$array[$j, $i];
    }
}