Search code examples
rsimulationmontecarlo

Monte Carlo Simulation with Replacement Based On Sum of A Column


I am trying to simulate an unlikely situation in a videogame using a Monte Carlo simulation. I'm extremely new at coding and thought this would be a fun situation to simulate.

There are 3 targets and they are being attacked 8 times independently. My problem comes with how to deal with the fact that one of the columns cannot be attacked more than 6 times, when there are 8 attacks.

I would like to take any attack aimed at column 2 select one of the other 2 columns at random to attack instead, but only if column 2 has been attacked 6 times already.

Here is my attempt to simulate with 5000 repeats, for example.

#determine number of repeats
trial <- 5000

#create matrix with a row for each trial
m <- matrix(0, nrow = trial, ncol = 3)

#The first for loop is for each row
#The second for loop runs each attack independently, sampling 1:3 at random, then adding one to that position of the row.
#The function that is called by ifelse() when m[trial, 2] > 6 = TRUE is the issue.

for (trial in 1:trial){
  for (attack in 1:8) {
    target <- sample(1:3, 1)
    m[trial, target] <- m[trial, target] + 1
    ifelse(m[trial, 2] > 6, #determines if the value of column 2 is greater than 6 after each attack
           function(m){
             m[trial, 2] <- m[trial, 2] - 1 #subtract the value from the second column to return it to 6
             newtarget <- sample(c(1,3), 1) #select either column 1 or 3 as a new target at random
             m[trial, newtarget] <- m[trial, newtarget] + 1 #add 1 to indicate the new target has been selected
             m}, #return the matrix after modification
           m) #do nothing if the value of the second column is <= 6
  }
}

For example, if I have the matrix below:

> matrix(c(2,1,5,7,1,0), nrow = 2, ncol = 3)
     [,1] [,2] [,3]
[1,]    2    5    1
[2,]    1    7    0

I would like the function to look at the 2nd line of the matrix, subtract 1 from 7, and then add 1 to either column 1 or 3 to create c(2,6,0) or c(1,6,1). I would like to learn how to do this within the loop, but it could be done afterwards as well.

I think I am making serious, fundamental error with how to use function(x) or ifelse.

Thank you.


Solution

  • Here's an improved version of your code:

    set.seed(1)
    
    trial <- 5000
    
    #create matrix with a row for each trial
    m <- matrix(0, nrow = trial, ncol = 3)
    
    #The first for loop is for each row
    #The second for loop runs each attack independently, sampling 1:3 at random, then adding one to that position of the row.
    #The function that is called by ifelse() when m[trial, 2] > 6 = TRUE is the issue.
    
    for (i in 1:trial){
        for (attack in 1:8) {
            target <- sample(1:3, 1)
            m[i, target] <- m[i, target] + 1 
            #determines if the value of column 2 is greater than 6 after each attack
            if(m[i, 2] > 6){ 
                #subtract the value from the second column to return it to 6
                m[i, 2] <- m[i, 2] - 1 
                #select either column 1 or 3 as a new target at random
                newtarget <- sample(c(1,3), 1)
                #add 1 to indicate the new target has been selected
                m[i, newtarget] <- m[i, newtarget] + 1 
                }   
            }   
        }   
    
    # Notice the largest value in column 2 is no greater than 6.
    apply(m, 2, max)
    

    set.seed is used to make the results reproducible (usually just used for testing). The ifelse function has a different purpose than the normal if-else control flow. Here's an example:

    x = runif(100)
    ifelse(x < 0.5, 0, x)
    

    You'll notice any element in x that is less than 0.5 is now zero. I changed your code to have an if block. Notice that m[i, 2] > 6 returns a single TRUE or FALSE whereas in the small example above, x < 0.5 a vector of logicals is returned. So ifelse can take a vector of logicals, but the if block requires there be only a single logical.

    You were on the right track with using function, but it just isn't necessary in this case. Often, but not always, you'll define a function like this:

    f = function(x)
        x^2
    

    But just returning the value doesn't mean what you want is changed:

    x = 5
    f(5) # 25
    x    # still 5
    

    For more on this, look up function scope in R.

    Lastly, I changed the loop to be i in 1:trial instead of trial in 1:trial. You probably wouldn't notice any issues in your case, but it is better practice to use a separate variable than that which makes up the range of the loop.

    Hope this helps.

    P.S. R isn't really known for it's speed when looping. If you want to make things goes faster, you'll typically need to vectorize your code.