# Linearity of Expectations - Repeated Card Draw

Scenario: You are given a standard deck of cards and told to draw one card at a time from the deck and record the color of that card. Continue until all the cards have been drawn without replacing any cards. Let the random variable $Z$ represent the number of times that you see a red card followed by a black card.

Question: What is $\mathbb{E}(Z)$?

Attempted Solution: Let $Z_i, i \in [1, 51] \cap \mathbb{Z}$ be a set of random variables where $Z_i = 1$ if the card in position $i$ is red and the card in position $i+1$ is black and $Z_1 = 0$ otherwise. Since $Z = \sum_{i = 1}^{51}Z_i$, we have the following by linearity of expectations:

$$\mathbb{E}(Z) = \sum_{i = 1}^{51}\mathbb{E}(Z_i)$$

Moreover, since the probability of drawing a red card followed by a black card without replacement from a standard deck is $\frac{26}{52} \cdot \frac{26}{51}$, we have: $$\mathbb{E}(Z) = \sum_{i = 1}^{51}\bigg(\frac{26}{52} \cdot \frac{26}{51}\bigg) = 51 \cdot \bigg(\frac{26}{52} \cdot \frac{26}{51}\bigg) = \frac{26^2}{52} = 13$$

Issue: The fact that linearity of expectations holds even when the random variables are dependent seemed like witchcraft to me, so I wrote a Python script to determine an experimental value of Z. This program gives me that on average $Z \approx 11.8$. Am I misinterpreting how to apply the linearity of expectations in this problem?

Here's the script I used:

from random import *

output_values = [];

n = 0;

while n < 10000:

#draw cards
red = 26;
black = 26;

draw = [None] * 52;

i = 0;

while i < len(draw):
if red > 0 and black > 0:
choose = randint(0,1);
if choose == 0:
draw[i] = "R";
red = red - 1;
else:
draw[i] = "B"
black = black - 1;
elif red == 0:
draw[i] = "B"
black = black - 1;
elif black == 0:
draw[i] = "R"
red = red - 1;
i += 1;

#count appearances of red followed by black
count_rb = 0;

j = 0

while j < (len(draw) - 1):
if draw[j] == "R" and draw[j+1] == "B":
count_rb += 1;
j += 1;

#add number of apperances to output_values array
output_values.append(count_rb);

n += 1;

average_value = (sum(output_values))/len(output_values);
print("Experimental Value:", average_value);

expected_value = 51 * (26*26)/(52*51);
print("Expected Value: ", expected_value);

• Could you add the Python script into the question (maybe the error is there) – Juho Kokkala Oct 16 '17 at 6:28
• I just added the script to the body of the question – ben-fogarty Oct 16 '17 at 14:24
• How does your script model drawing without replacement from the deck of cards? It doesn't appear to use the correct probabilities. To understand the problem, consider how your code works with four cards. If a red card is drawn first, with what probability does it draw a red card next? What should the correct probability be? – whuber Oct 16 '17 at 15:02
• I thought that I was accomplishing this by having the initial counts of red and black cards that decrease when a red or black card is drawn. Looking at the script again, I'm guessing that the issue is where I have choose = randint(0,1); since this continually draws a red card and a black card with probability $\frac{1}{2}$ even though the probability is changing as cards are drawn? – ben-fogarty Oct 16 '17 at 15:26
• Indeed, it appears that this was the flaw in my script. I've posted a corrected script as an answer below. – ben-fogarty Oct 16 '17 at 15:50

As suggested by @whuber and @Juho Kokkala, the issue was with my simulation and not with my calculation. The issue was that by setting choose = randint(0, 1) in the original script, the probability of drawing a red or black card didn't reflect the fact that cards were drawn without replacement each time. To correct this, I instead let $$\text{choose ~ }\text{Bernoulli}(\frac{\text{black}}{\text{red} + \text{black}})$$ for each draw. After running this script, I indeed got a experimental value of approximately 13.

For the curious, I present my entire revised script below. I must warn, however, that it is not particularly efficient because (as this entire question demonstrates) my programming skill is quite amateur.

from scipy.stats import bernoulli

output_values = []

n = 0

while n < 5000:
red = 26
black = 26

draw = [None] * 52

i = 0

while i < len(draw):
if red > 0 and black > 0:
choose = bernoulli.rvs(black/(red + black))
if choose == 0:
draw[i] = "R"
red = red - 1
else:
draw[i] = "B"
black = black - 1
elif red == 0:
draw[i] = "B"
black = black - 1
elif black == 0:
draw[i] = "R"
red = red - 1
i += 1

count_rb = 0

j = 0

while j < (len(draw) - 1):
if draw[j] == "R" and draw[j+1] == "B":
count_rb += 1
j += 1
output_values.append(count_rb)

n += 1

average_value = sum(output_values)/len(output_values)
print("Experimental Value:", average_value)

expected_value = 51 * (26*26)/(52*51)
print("Expected Value: ", expected_value)

• (+1) This code for an R simulation might suggest additional ways to improve the calculation: mean(replicate(5e3,(function(x)sum(!x[-1]&x[-length(x)]))(sample(rep(0:1,26))))) – whuber Oct 16 '17 at 17:14