Expected value of sum of cards

Question

If each card on a regular 52 deck card has points that corresponds to their number (like 2 of hearts is 2 points, 7 of clubs is 7 points), the Jack, Queen, King each being 10 points and you keep drawing without replacement until sum of all points is 10 or more....what's the mean of sum of points?

Answer 1

By your definition, you have $16$ cards ($10$, $\text{J}$, $\text{K}$, $\text{Q}$) that are worth $10$ points, so with probability $16/52$ you get $10$ points in a single draw. Since $9+\text{anything}=10$, then if we take into consideration that there is $4$ nines, than we instantly know that with probability greater than $20/52$ you finish in two draws. However, return of two draws is simple to obtain by enumerating all $52 \choose 2$ combinations of card pairs and summing their scores.

unique_cards <- c(1:10, 10, 10, 10)  # A, 1, 2, ..., 10, J, Q, K
unique_cards <- rep(unique_cards, 4) # each appears 4 times

comb <- combn(unique_cards, 2) # take all possible combinations of card pairs

what gives $79\%$ probability of obtaining score of at least $10$ in two draws

> sum(colSums(comb) >= 10)/choose(52, 2) # accepted / all combinations
[1] 0.7888386

Lazy solution for more than two draws can be obtained by a simple simulation, where whole deck is shuffled and then cards are drawn until their total score is at least $10$.

set.seed(123)

sim <- function(target = 10) {
  res <- cumsum(sample(unique_cards)) # shuffle, draw and sum
  n <- which.max(res >= target)       # take first score >= 10
  c(sum = res[n], n = n)
}

R <- 1e4
res <- replicate(R, sim())

and the result is that on average you have to draw two cards and the average total score is $12.77$

> apply(res, 1, summary)
          sum     n
Min.    10.00 1.000
1st Qu. 10.00 1.000
Median  12.00 2.000
Mean    12.77 1.946
3rd Qu. 15.00 2.000
Max.    19.00 6.000

Moreover, as expected, with approximately $30\%$ probability you finish with one draw, but with $79\%$ probability you finish with two draws and you rarely get over three draws:

 > cumsum(table(res[2,])/R)
     1      2      3      4      5      6 
0.3006 0.7910 0.9653 0.9968 0.9999 1.0000 

Answer 2

The question asks for the "mean of the sum of points." Because the target is 10 and no value exceeds 10, this is the mean of a distribution defined on the ten integers $10, 11, \ldots, 10+10-1$.

It takes about the same amount of computing power to work out this distribution, exactly, as it does to perform a small simulation. Here's how.

The deck is determined by the numbers of cards of each point. Let there be $k_1 = 4$ cards worth one point (the aces), $k_2$ worth two points, and so on. Writing $n=10$, the vector $\mathbf{k}=(k_1,k_2, \ldots, k_n)$ for a standard deck is

$$\mathbf{k} = (4,4,4,4,4,4,4,4,4,16).$$

When we draw a card from this deck, we remove a card worth $i$ points with probability

$$P_{\mathbf{k}}(i) = \frac{k_i} {k_1+k_2+\cdots+k_n}.$$

The deck changes afterwards: $k_i$ is reduced to $k_i-1$. Let's indicate the new vector with the notation

$$S_i(\mathbf{k}) = (k_1, k_2, \ldots, k_{i-1}, k_i - 1, k_{i+1}, \ldots, k_n).$$

Let $f_{\mathbf{k}}(t, s)$ be the distribution of the sum of points when the target is $t$ starting with a sum of $s$ points. We wish to find $f_{\mathbf{k}}(10, 0)$.

The possible draws are described by $n=10$ distinct, non-overlapping events: event $i$ consists of drawing a card worth $i$ points. When that happens, any starting sum $s$ is increased to $s+i$, the target we need to reach is reduced to $t-i$, and the deck is changed to $S_i(\mathbf{k})$. The Law of Total Probability tells us to sum the chances over all these events. Thus,

$$f_{\mathbf{k}}(t, s) = \sum_{i=1}^n P_{\mathbf{k}}(i)f_{S_i(\mathbf{k})}(t-i, s+i).\tag{*}$$

Certainly when $s$ exceeds the original target ($10$) there's nothing left to figure out: the drawing terminates and the distribution will be 100% on the total value $s$. Equivalently, when the target is $0$ or negative then we should put all the probability on whatever value $s$ currently has, because the target obviously has been met.

These considerations give an effective recursive formula $(*)$ for $f$. Because $t$ decreases by at least $1$ with each iteration, it is guaranteed to terminate within $10$ draws. (It will actually terminate within $7$ draws, because soon all the aces would be used up.) This is short. (Just 4389 calls to $f$ are needed and only 446 of those have to perform the summation in $(*)$.) With double-precision computations the answer takes only a tenth of a second to obtain. It gives these probabilities for the final values $10, 11, \ldots, 19$ (which have been rounded for ease of reading):

0.37906 0.09108 0.08503 0.08203 0.07512 0.07132 0.06356 0.05877 0.04995 0.04408

Their expectation is $$10\times 0.37906 + 11\times 0.09108 + \cdots + 19\times 0.04408 = 12.7534.$$

---

An R implementation of $f$ is given below. First, though, its results can be checked by simulation. The game is played $10^4$ times, the final sums are tallied, and those tallies are compared to the preceding results with a $\chi^2$ test. (It is applicable and accurate because the smallest expected value of any cell count is a very large $440.8$.)

set.seed(17)
deck <- c(rep(4, 9), 16) # deck[i] counts cards of value `i`.
deck.long <- unlist(sapply(1:length(deck), function(i) rep(i, deck[i])))
sim <- replicate(1e4, {
  x <- cumsum(sample(deck.long, 10))
  x[which(x >= 10)[1]]
})
y <- table(sim)
z <- c(rep(0, 10), y/sum(y))
rbind(x, z)
chisq.test(y, p=x[-(1:10)])

The output is

    Chi-squared test for given probabilities

data:  y
X-squared = 3.8856, df = 9, p-value = 0.9188

The large p-value demonstrates consistency between the simulation and the theoretical answers.

Here is how $f$ can be computed.

f <- function(deck, total=0, target=10, maximum=20) {
  x <- rep(0, maximum)
  if (target <= 0) {
    x[total+1] <- 1
    return(x)
  }
  n <- sum(deck)
  x <- sapply(1:length(deck), function(i) {
    k <- deck[i]
    if (k <= 0) return (x)
    d <- deck
    d[i] <- d[i] - 1
    k/n * f(d, total + i, target - i, maximum)
  })
  return(rowSums(x))
}
x <- f(deck)
round(x[-(1:10)], 5)
sum(x * (1:length(x)-1)) # Expected value

---

Incidentally, the same techniques--with only the tiniest modifications (which I leave to interested readers)--will answer the question "what is the distribution of the number of draws in the game?" The answer (again in double precision) is the vector of probabilities corresponding to 1, 2, ..., 7 draws:

3.076923e-01 4.811463e-01 1.762293e-01 3.176286e-02 3.029212e-03 1.381895e-04 1.866540e-06

A similar simulation--this time of a million iterations, because these probabilities get so small--produces these observed frequencies:

     1      2      3      4      5      6      7 
306897 481652 176242  32052   3019    135      3 

These do not differ significantly from the computed values.

Answer 3

I ran the simulation is C# .NET
I am getting 12.75x consistently

private static decimal AvgSumToTen()
{
    Int32 loops = Int32.MaxValue / 10;
    //loops = 100;
    Random rand = new Random();
    int thisSum;
    ulong ttl = 0;
    int b;
    int bTenRaw;
    int bTen;
    HashSet <int> values = new HashSet<int>();
    for(Int32 i = 0; i < loops; i++)
    {
        thisSum = 0;
        values.Clear();
        while (thisSum < 10)
        {
            b = rand.Next(0, 52);
            if (values.Contains(b))
                continue;
            values.Add(b);
            bTenRaw = b % 13 + 1;
            bTen = (bTenRaw >= 10) ? 10 : bTenRaw;
            //Debug.WriteLine("bTen "+ bTen);
            thisSum += bTen;
        }
        //Debug.WriteLine("thisSum " + thisSum + Environment.NewLine);
        ttl += (ulong)thisSum;
        if (ttl > (ulong.MaxValue - 100))
            Debug.WriteLine("ttl > (ulong.MaxValue - 100)" + thisSum);
    }
    decimal answer = (decimal)ttl / (decimal)loops;
    Debug.WriteLine("answer " + answer.ToString("N4"));
    return answer;
}