simulation of a vector at different points using a for loop in R

Question

I am having some problems trying to write a set of simulations in which the starting value of my variable is a vector to be evaluated at different times instead of a single observation. One could run the simulations one by one by changing the starting point and then regroup all the results into another matrix; however, I am dealing approximately with 2000 cases, thus a more elegant and faster way is needed.

This is not the real code I am using but an example of the same problem, let's say:

time<-c(65,130,195,260)  #in days

simulation<-matrix(numeric(10*4),10,4)  

  #vessel matrix containing the desired number of simulations for 1 initial case

results<-NULL

initial<-.5   

   #This would be one case and I've got a vector with approx.2000

   #loop

for(i in 1:10){
  simulation[i,] = (initial*exp(-a*(time/260)) + 
                   b*(1-exp(-a*(time/260)))    +
                   c*(j/260))*rnorm(1)
}

results<-colMeans(simulation)


  # With this I end up with a row vector of 4 entries containing the average of 
  # the simulations for the first case at four dates.

how would you replace initial with a vector, let's say initial=seq(from=0, to=10, by=.5), where I could end up with a matrix of 20X4 where each row still contains the average of 10 simulations for each date?

Answer 1

So, you can use replicate instead of a for loop, wrap it in a function, and then you're done.

time<-c(65,130,195,260) 
initial<-.5   

# I assume a,b,c,j are constants
a<-b<-c<-j<-1

# Use replicate instead of a four loop to generate your matrix directly:
mat<-replicate(10,(initial*exp(-a*(time/260)) + b*(1-exp(-a*(time/260))) + c*(j/260))*rnorm(1))
# Get the means for all four variables (note matrix is transposed)
rowMeans(mat)

# Wrap this in a function:
f<-function(initial) {
  mat<-replicate(10,(initial*exp(-a*(time/260)) + b*(1-exp(-a*(time/260))) + c*(j/260))*rnorm(1))
  rowMeans(mat)
}

# Now repeat for a variety of "initial"
sapply(seq(from=0, to=10, by=.5),f)

# [,1]       [,2]      [,3]      [,4]       [,5]      [,6]        [,7]       [,8]
# [1,] 0.04131179 -0.1703167 0.2165282 0.4518704 -0.4305232 1.0390553 -0.09262142 -0.3767832
# [2,] 0.07293559 -0.1941923 0.2165282 0.4239343 -0.3889186 0.9154405 -0.08016291 -0.3217918
# [3,] 0.09756422 -0.2127866 0.2165282 0.4021776 -0.3565170 0.8191691 -0.07046022 -0.2789645
# [4,] 0.11674502 -0.2272679 0.2165282 0.3852335 -0.3312825 0.7441930 -0.06290376 -0.2456105

Answer 2

time<-c(65,130,195,260)  #in days

simulation<-matrix(numeric(10*4),10,4)  
#vessel matrix containing the desired number of simulations for 1 initial case

results<-NULL

#initial<-.5   

initial=seq(from=0.5, to=10, by=.5) # Changed This to be vector


#This would be one case and I've got a vector with approx.2000

#loop
## Writing to a function to help in reusing the code
MM <- function(initial){

for(i in 1:10){
  simulation[i,] = (initial*exp(-10*(time/260)) + 
                      20*(1-exp(-10*(time/260)))    +
                      20*(1000/260))*rnorm(1)
}

results<-colMeans(simulation)
return(results)
}

sapply(initial, FUN=MM) # This prints the colmeans for each inputted element in the vector initial

#Output:
#                 [,1]            [,2]            [,3]           [,4]
#[1,] 0.04560954756720 -7.358191844904 -8.098682864513 3.665639204753
#[2,] 0.04631255769454 -7.468652842648 -8.217009769192 3.717726620815
#[3,] 0.04637026427987 -7.477720033496 -8.226722633000 3.722002216291
#[4,] 0.04637500112485 -7.478464313845 -8.227519913413 3.722353178540
#               [,5]           [,6]           [,7]           [,8]
#[1,] 18.98057953410 21.03382334801 11.96959043580 63.54595548535
#[2,] 19.24268217595 21.31585924907 12.12529951449 64.34721302134
#[3,] 19.26419687094 21.33901016563 12.13808089400 64.41298424507
#[4,] 19.26596290465 21.34091050858 12.13913005352 64.41838307588
#                [,9]          [,10]          [,11]           [,12]
#[1,] -20.01917307564 12.56144497430 27.31607000721 -17.70877647896
#[2,] -20.26360373127 12.70980687528 27.62780928081 -17.90382189201
#[3,] -20.28366782130 12.72198516172 27.65339839866 -17.91983219447
#[4,] -20.28531478210 12.72298481635 27.65549888136 -17.92114640013
#               [,13]           [,14]          [,15]           [,16]
#[1,] -1.163127428493 -13.07282052205 20.91363938663 -15.43763356270
#[2,] -1.175475358083 -13.20640615074 21.11903955938 -15.58312496123
#[3,] -1.176488937866 -13.21737152689 21.13589983227 -15.59506762248
#[4,] -1.176572137562 -13.21827161978 21.13728380775 -15.59604793581
#              [,17]           [,18]           [,19]           [,20]
#[1,] 4.129982905656 -32.29365824471 -26.80599719664 -13.92279757936
#[2,] 4.167268037791 -32.57240758567 -27.02676696018 -14.03195652202
#[3,] 4.170328587811 -32.59528872494 -27.04488884592 -14.04091683368
#[4,] 4.170579813055 -32.59716692322 -27.04637638088 -14.04165234085

You should use apply family of functions if you want to use a code over a vector or list. Hope this helps !