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);

 A: 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)

