0

I'm trying to simulate 1000 trials of the following situation to get a mean output of cost. In addition to what's included below, each simulation would need to stop if the totalTime >= 6. I struggle a lot with for loops and simulations so I would appreciate it if anyone can give me any tips for creating this simulation. Thanks!

cost <- 0
devices <- 0
totalTime <- 0
unitCost <- 100

failureTime <- rgamma(1, shape=2, scale=.5)
if (failureTime>1) {
  cost <- cost + unitCost
  totalTime <- totalTime + failureTime
} else {
  cost <- cost
  totalTime <- totalTime + failureTime
}
  • 2
    Your notion of "stop" seems to suggest that in addition to repeating this full-simulation 1000 times, you want each simulation to loop until a condition is met. If you have a basic 1-full-simulation working somewhere, then you can use `replicate(1000, { your_code_here; })`. – r2evans Jun 18 '21 at 20:08

2 Answers2

1

You can use replicate like below

replicate(1000, {
    cost <- 0
    totalTime <- 0
    unitCost <- 100
    while (totalTime < 6) {
        failureTime <- rgamma(1, shape = 2, scale = .5)
        cost <- cost + unitCost * (failureTime > 1)
        totalTime <- totalTime + failureTime
    }
    cost
})

Note that you need a while loop to judge if totalTime reaches 6, and cost (or totalTime as well) is output once the while loop is terminated.

If you want data.frame output, you can try the code below

out <- do.call(
    rbind,
    replicate(1000,
        {
            cost <- 0
            totalTime <- 0
            unitCost <- 100
            while (totalTime < 6) {
                failureTime <- rgamma(1, shape = 2, scale = .5)
                cost <- cost + unitCost * (failureTime > 1)
                totalTime <- totalTime + failureTime
            }
            data.frame(cost = cost, totalTime = totalTime)
        },
        simplify = FALSE
    )
)

which gives

> head(out)
  cost totalTime
1  200  6.117942
2  200  6.010357
3  400  6.255933
4  400  6.082638
5  300  7.170587
6  200  7.046343
ThomasIsCoding
  • 96,636
  • 9
  • 24
  • 81
1

I think it's important to first have a full-running simulation work 1 time. To do that and honor your "stop if the totalTime >= 6", you need something like a while loop.

I suggest writing a function that does this one thing very well and returns the two variables you want (cost and total time). Once you have that, it's trivial to replicate it.

simfunc <- function(shape = 2, scale = 0.5, unitCost = 100) {
  cost <- 0
  totalTime <- 0
  while (totalTime < 6) {
    failureTime <- rgamma(1, shape = shape, scale = scale)
    if (failureTime > 1) {
      cost <- cost + unitCost
    }
    totalTime <- totalTime + failureTime
  }
  c(cost = cost, totalTime = totalTime)
}

set.seed(42)
simfunc()
#      cost totalTime 
#    200.00      6.86 
simfunc()
#      cost totalTime 
#    400.00      6.79 
simfunc()
#      cost totalTime 
#    300.00      6.95

From here, we can replicate it easily:

set.seed(42)
replicate(10, simfunc())
#             [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7] [,8]   [,9]  [,10]
# cost      200.00 400.00 300.00 300.00 200.00 300.00 300.00  200 400.00 300.00
# totalTime   6.86   6.79   6.95   6.62   6.34   6.81   7.14    6   6.73   6.75

(To summarize, store that in a variable and aggregate the row you need. You might find it easier to visualize if you transpose it to be columnar, and perhaps convert to a data.frame, but that's only for your aesthetic benefit, the stats will be the same regardless.)

Also, I'm inferring that you may want totalTime, since often simulations like this not only want to know cost over the system's life-cycle, but also the mean-time between fails and some analysis of that. In this case, it appears you've already figured that out (based on your scale= and shape= parameters), but it can be interesting anyway.

r2evans
  • 141,215
  • 6
  • 77
  • 149