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)

 A: It's not clear from your question whether you are sampling with or without replacement.  In the former case the probability of winning is $p_*=0.52$ and in the latter case the probability of winning is $p_{**}=0.40$.  (As an exercise, try to derive these probabilities from your problem.)  Consequently, the distribution of the profit and its mean and variance under either scenario is:
$$\begin{matrix}
\ \ \pi_* \sim 100 \cdot \text{Bin}(17, 0.52) - 80*17
& & & \ \mathbb{E}(\pi_*) = -476
& & & \ \mathbb{V}(\pi_*) = 42432, \\[6pt]
\pi_{**} \sim 100 \cdot \text{Bin}(17, 0.40) - 80*17
& & & \mathbb{E}(\pi_{**}) = -680
& & & \mathbb{V}(\pi_{**}) = 40800. \\[6pt]
\end{matrix}$$
Both of these results can easily be confirmed by simulation.  In R we can use the sample.int function to directly sample from the set under consideration.  Here I will construct a function to simulate the profit for this kind of problem.  For generality I will take the specific values in your game and set them as default inputs that allow alternative specification in the function.
#Create function to simulate the profit
simulate.profit <- function(sims, games = 17, 
                            objects = 5, draws = 2, replace = FALSE,
                            cost = 80, winnings = 100) {
  
  #Determine outcomes of games and profit
  WIN    <- matrix(FALSE, nrow = sims, ncol = games)
  PROFIT <- rep(0, sims)
  for (i in 1:sims) {
    for (j in 1:games) {
      SAMPLE   <- sample.int(objects , size = draws, replace = replace)
      WIN[i,j] <- (sum(SAMPLE) %% 2 == 0) } 
    PROFIT[i] <- winnings*sum(WIN[i,]) - games*cost }
  
  #Return profit simulations
  PROFIT }

Using this function we can easily simulate the profit values under either sampling-with-replacement or sampling-without-replacement.  The simulations below confirm the means and variances calculated above.
#Simulate using sampling-without-replacement
set.seed(1)
SIMS <- 10^7
PROFIT.SIMS1 <- simulate.profit(sims = SIMS, replace = TRUE)

#Compute the mean
mean(PROFIT.SIMS1)
[1] -476.0598

#Compute the variance
var(PROFIT.SIMS1)
[1] 42432.42

###########################################################
#Simulate using sampling-with-replacement
set.seed(1)
SIMS <- 10^7
PROFIT.SIMS2 <- simulate.profit(sims = SIMS)

#Compute the mean
mean(PROFIT.SIMS2)
[1] -680.0054

#Compute the variance
var(PROFIT.SIMS2)
[1] 40802

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