Probability game - Simulation don't match the manual solution Two numbers are bring chosen by random from the set {1,2,3,4,5}. If the sum of the two numbers is even, you win 100 dollars, otherwise you win nothing. In order to participate in the game, you pay $80. What is the expected value and variance of the profit after 17 games ?
I solved this one analytically and with Monte Carlo simulation and got different results, I wonder where my mistake it.
My Solution:
The probability of success in one round is:
$p=\frac{\binom{2}{2}\binom{3}{0}+\binom{2}{0}\binom{3}{2}}{\binom{5}{2}}=0.6$
The number of successes is:
$X:Bin(17,0.6)$
Therefore:
$E(X)=10.2 , V(X)=4.08$
The profit is defined as:
$Y=100X-1360$
And thus:
$E(Y)=-340 , V(Y)=40800$
My R code is below, and it yields different results , e.g.:
$E(Y)=-480.38 , V(Y)=42751.21$
The differences are too big in my opinion, so there must be a mistake somewhere.
Can you assist me identifying the mistake ?
n = 5000
successes = rep(0,n)
for (j in 1:n)
{
  result = rep(0,17)
  for (i in 1:17)
  {
    game.result = sample(1:5,2,replace = T)
    if (sum(game.result)%%2==0)
    {
      result = 1
    }
  }
  successes[j] = sum(result)
}

successes
mean(successes)
var(successes)

profit = 100*successes-1360
mean(profit)
var(profit)


As others noticed, you made few mistakes in your assumptions and the code. In your simulation you sample with replacement. The description of the problem says nothing about drawing without replacement, so it's up to you how you read it. In such a case, the probability $p$ is
# list all the possibilities
outcomes <- matrix(NA, 5, 5)
for (i in 1:5) {
    for (j in 1:5) {
        outcomes[i,j] <- i + j
    }
}
# how many there are?
m <- prod(dim(outcomes))
# how many are winning?
k <- sum(outcomes %% 2 == 0)
p <- k / m
## 13 / 25 = 0.52

n <- 17
(n * p * 100) - (17 * 80)
## -476

Also there are bugs in your code. When writing a simulation code, it is usually a good idea to make it as simple as possible, so that it literally matches the problem. In such a case, you have the lowest risk of incorrect results due to the bug.
sim <- function() {
    profit <- 0
    for (i in 1:17) {
        outcome <- sample(1:5, 2, replace = T)
        profit <- profit - 80
        if (sum(outcome) %% 2 == 0) {
            profit <- profit + 100
        }
    }
    return(profit)
}

result <- replicate(50000, sim())

summary(result)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1360.0  -660.0  -460.0  -475.2  -360.0   340.0 
var(result)
## 42469.426952539

of using more high-level simulation:
result <- (rbinom(50000, size=n, prob=p) * 100) - (17 * 80)

summary(result)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1260.0  -660.0  -460.0  -476.4  -360.0   340.0 
var(result)
## 42303.3516630333

As you can see, they match.
