saving output from a nested for-loop in R

Question

I'm trying to figure out how to store the output of a simulation that I will run for many iterations. I have a simplified example here where I'm simulating changes in winter temperature over a 50-year period, and I'd like to do this 3 times:

# function to draw the winter temperature in a particular year
draw_winter_temperature <- function(){
  intercept <- winterIntercept
  beta1 <- rnorm(1, winterBeta, winterBetaSD)
  val = rnorm(1, intercept + beta1 * years[i], winterSigma)
  return(val)
}

# parameters for winter temperature function
winterIntercept <- -0.633
winterBeta <- 0.08
winterBetaSD <- 0.04
winterSigma <- 1.30 

# number of iterations
iter = seq(1,3,1)

# number of years to run in each simulation
years <- seq(1,50,1)

# matrix to store 50 winter temperatures for each iteration? Not sure best approach
temps <- matrix(NA, nrow = length(years), ncol = length(iter))

# start of nested for loop - loop over each iteration
for (i in seq_along(iter)) {

  # within each iteration, loop over each year
  for (j in seq_along(years)) {

  # draw winter temp. anomaly
  winterTemp <- draw_winter_temperature()

  # store winter temps, not sure what to do here
  temps <- winterTemp

  }
}

I guess I'm a little unsure about the indexing here to save the output in my empty matrix (or to save it in another format). I should say that there are other functions and data being stored within my loop over years, so I'm generally confused about the best way to store data from each of the iterations.

Answer 1

I think you need the [j, i] subset in temps. Also, redefine the function to move across the years argument:

> rm(list = ls())
> 
> # parameters for winter temperature function
> winterIntercept <- -0.633
> winterBeta <- 0.08
> winterBetaSD <- 0.04
> winterSigma <- 1.30 
> 
> # number of iterations
> iter = seq(1,3,1)
> 
> # number of years to run in each simulation
> years <- seq(1,50,1)
> 
> # matrix to store 50 winter temperatures for each iteration? Not sure best approach
> temps <- matrix(NA, nrow = length(years), ncol = length(iter))
> 
> draw_winter_temperature <- function(i) {
+   intercept <- winterIntercept
+   beta1 <- rnorm(1, winterBeta, winterBetaSD)
+   val = rnorm(1, intercept + beta1 * years[i], winterSigma)
+   return(val)
+ }
> 
> 
> # start of nested for loop - loop over each iteration
> for (i in seq_along(iter)) {
+   
+   # within each iteration, loop over each year
+   for (j in seq_along(years)) {
+     
+     # draw winter temp. anomaly
+     winterTemp <- draw_winter_temperature(j)
+     
+     # store winter temps, not sure what to do here
+     temps[j, i] <- winterTemp
+     
+   }
+ }
> 
> temps
             [,1]        [,2]        [,3]
 [1,]  0.07484416 -2.09659305  1.36146643
 [2,]  0.91077388  0.91630545  0.15919461
 [3,]  2.07646459  1.47967364 -0.80436293
 [4,] -1.06727442 -1.35760417 -2.23396119
 [5,]  0.67833459 -0.49733006 -0.08458830
 [6,] -2.94754201 -0.83298461 -0.09547269
 [7,] -0.27718212 -1.26709144  1.07166328
 [8,] -2.48913267  1.15023058 -0.45768205
 [9,]  0.81201831  1.73791432 -1.12876304
[10,] -1.26716059 -0.59580422 -1.28796699
[11,]  1.73244202  1.45400888 -0.13374732
[12,]  0.06181381  1.42056333 -0.28493854
[13,]  1.32197105 -0.26822322 -0.64071024
[14,]  0.31681581  2.17959219  1.05139421
[15,] -0.40649165  0.40584711 -1.24801466
[16,] -0.37852770 -0.01912727 -1.57738425
[17,]  1.07342597 -0.37121912 -1.78913387
[18,]  2.37652192  2.34135284  1.57626300
[19,]  2.26526300 -0.22457651  1.22608845
[20,]  0.43476700  1.20513841  3.15556064
[21,]  2.59591314 -0.30772162 -0.18688845
[22,]  1.06465892  1.71525510 -1.96024033
[23,]  0.33484399  1.81225896  2.11832405
[24,]  1.90063594  3.51556209  1.07306117
[25,] -0.01940422  2.28460208  3.21963212
[26,]  3.28816991  2.95440512  3.24107428
[27,]  0.96753190 -0.10513793 -1.21244216
[28,]  2.54233795 -0.71844973  1.81812858
[29,]  2.17954621  2.98536138  3.66748584
[30,]  0.81877427 -1.35965678  2.56338093
[31,]  0.29897273  1.11832993  0.87086081
[32,]  3.41605791  4.72562882  4.07046433
[33,]  3.48441756 -0.09328612  1.64075691
[34,] -1.45902220 -0.38229884  4.62577573
[35,]  0.76568048  3.05618914 -2.04568644
[36,]  0.51939827  0.46442133  2.45079102
[37,]  3.31122478  3.72294771  3.79475270
[38,]  3.22139644  2.33785724 -0.38787472
[39,]  2.95662591  1.14475461  2.97768620
[40,]  0.41263989  3.12365715 -0.05141565
[41,]  1.38138984  3.76526281  1.96114126
[42,] -0.18932163  3.05150559  1.01252868
[43,]  2.60596068  4.58785954 -0.62515705
[44,]  4.11088547  1.61733269  2.51478319
[45,]  2.54634136  5.03258017  3.69739035
[46,]  3.90668226  5.10453707  9.44805817
[47,] -0.01088235 -0.72829163  5.23958084
[48,]  5.69817388  0.12111318  4.43908874
[49,]  0.62863949  1.44337538  6.37749772
[50,]  6.03095905 -1.06661804  2.96625244

It can also be done with nested sapply:

> sapply(seq(1, 3, 1), function(x){
+   sapply(seq(1, 50, 1), function(i) {draw_winter_temperature(i)})
+ })
             [,1]       [,2]        [,3]
 [1,] -0.89724777 -0.3248484 -0.90899374
 [2,] -1.94179327 -1.1282839  0.31149541
 [3,] -0.42025521  0.7724711  0.74888190
 [4,] -1.33109020  1.9129040  0.87848969
 [5,] -0.47958978  0.0847167  0.41101536
 [6,] -0.79725942 -0.3091579 -1.24232376
 [7,]  1.77518276 -0.1922316 -1.47421230
 [8,]  0.69493304 -0.2019549  0.13570145
 [9,] -1.30457205  0.3104541  0.42987422
[10,]  0.17036195 -1.7092747 -0.83657188
[11,]  0.32147101 -0.9007392 -1.19378838
[12,]  0.67253811  1.3546455  0.65988074
[13,]  1.50350007  0.5492711  1.79631547
[14,] -2.51666246 -0.9415391 -0.29810167
[15,] -0.40993616 -0.4815619 -0.64542730
[16,]  1.77118058  0.1315932 -0.75448968
[17,]  0.50048465  2.7956318 -0.87173764
[18,] -3.21569892 -0.3323585 -1.23359412
[19,]  0.31825702  0.3584551  0.06093396
[20,]  3.89144244 -2.4149095  1.78682882
[21,] -0.43390072  0.7793848 -0.10369739
[22,] -2.86730666 -0.7670008  1.52314298
[23,] -1.79228394  0.9572841  0.48042120
[24,] -0.43718475 -0.1275890  1.21647106
[25,] -1.15994958  1.6406028  4.60036997
[26,]  0.01317481  1.4994584  2.66739115
[27,]  4.17871392  0.7907429  2.73681377
[28,]  4.28581141  1.6371291  0.68691671
[29,]  1.93210763  2.9950647  4.21374229
[30,]  2.58606149  1.0610567 -0.85244266
[31,]  1.80716913  4.5783347  2.52572309
[32,] -0.59753770  1.4036211  4.67023114
[33,] -0.62456359  3.8556583  1.45738330
[34,]  3.97709081  1.9414286  2.13243487
[35,] -0.44078059  0.1845302  2.93266443
[36,]  3.95186775  0.6522162  1.33296220
[37,]  2.44364830  2.8547030  1.20660793
[38,]  4.41777409  1.9200504  2.71936027
[39,]  0.28099830  2.7876076  3.92335543
[40,]  3.40916089  0.1663190  0.42249808
[41,]  4.14691320  2.2004945  2.03032191
[42,]  0.61603992  4.6339235  2.82358688
[43,]  3.15565113  4.2119661  5.04334175
[44,]  2.63957253  5.9546210  2.77191218
[45,] -1.50822755  0.9417640  6.71999760
[46,]  1.29101617  3.2551090  1.59456129
[47,]  5.98411236  2.0040345  4.03826529
[48,]  2.33718663 -3.1029647  2.84654953
[49,] -1.77860722  4.4929658  2.93310143
[50,]  7.12512178  5.3284199  2.97240287

Hope it helps.

Answer 2

I'd avoid the loops altogether. Create a data frame with the iterations and years:

data1 <- data.frame(iter = rep(1:3, each = 50), 
                    year = rep(1:50, 3))

Define the function a little differently, so as all the default parameters are supplied:

draw_winter_temperature <- function(intercept = -0.633, 
                                    winterbeta = 0.08, 
                                    winterbetasd = 0.04,
                                    wintersigma = 1.30,
                                    year = 1){
  beta1 <- rnorm(1, winterbeta, winterbetasd)
  val   <- rnorm(1, intercept + beta1 * year, wintersigma)
  val
}

Use purrr::map_dbl() and dplyr::mutate() to create a column with sampled values for each year:

library(dplyr)
library(purrr)

data1 <- data1 %>% 
  mutate(winterTemp  = map_dbl(year, ~draw_winter_temperature(year = .)))