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)


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

