Speed of Simulations

We will compare different methods of simulating 100,000 realizations of the sample mean from samples of size 9 from an exponential distribution with rate, λ = 1.

N <- 100000
n <- 9

For reference, the method suggested in Simulation studies using R is

set.seed(123)
system.time(resRepl <- replicate(N, mean(rexp(n))))

   user  system elapsed 
  3.630   0.000   3.671

We see that on the machine I am using it takes a few seconds to produce the result, which is

str(resRepl)

 num [1:100000] 0.704 0.584 1.505 1.083 1.489 ...

set.seed(123)
resLoop1 <- numeric(N)        # preallocate the result of correct size
system.time(for (i in 1:N) resLoop1[i] <- mean(rexp(n)))

>    user  system elapsed 
  3.880   0.000   3.904

all.equal(resRepl, resLoop1)

[1] TRUE

set.seed(123)
resLoop2 <- numeric()                   # empty numeric vector
system.time(for (i in 1:N) resLoop2[i] <- mean(rexp(n)))

>    user  system elapsed 
 42.480   0.050  43.018

all.equal(resLoop2, resRepl)

[1] TRUE

set.seed(123)
resLoop3 <- NULL                        # empty list
system.time(for (i in 1:N) resLoop3 <- c(resLoop3, mean(rexp(n))))

>    user  system elapsed 
 40.930   0.060  41.073

str(resLoop3)

 num [1:100000] 0.704 0.584 1.505 1.083 1.489 ...

At one time it was important to take into consideration the overhead of calls to R functions and avoid putting too many functions inside a loop. That is less of an urgent consideration on modern computers. Although the interpreter overhead is still present on modern computers it is not as big an issue.

Nevertheless it is good to see how such a simulation could be done by creating all the random numbers from the exponential distribution in one call and then reducing each subsample. We arrange the n*N values into an n× N matrix. (We could make it N× n but there is a slight advantage in having subsamples in columns rather than rows.)

set.seed(123)
str(MM <- matrix(rexp(n*N), nr = n))

 num [1:9, 1:100000] 0.8435 0.5766 1.3291 0.0316 0.0562 ...

set.seed(123)
system.time(resColMeans <- colMeans(matrix(rexp(n*N), nr = n)))

   user  system elapsed 
  0.150   0.000   0.143

all.equal(resColMeans, resRepl)

[1] TRUE

resApply <- apply(MM, 2, mean)
all.equal(resApply, resRepl)

[1] TRUE

set.seed(123)
system.time(resApply <- apply(matrix(rexp(n * N), nr = n), 2, mean))

   user  system elapsed 
  4.000   0.000   4.058

The apply function is one of a family of functions that apply other functions to the elements of some structure. The most general of these is lapply, which applies a function to the elements of a list (or any other vector structure, including numeric vectors).

A characteristic of lapply is that it always returns a list. The sapply function is like lapply except that it tries to "simplify" the list that is returned.

In our case we just want to repeat the operation of simulating n random values and determining their mean and do that N times so the function we want to apply doesn't use its argument. Nevertheless, we must give it an argument, even though we don't use it, so that the sapply function stays happy.

set.seed(123)
system.time(resSapply <- sapply(1:N, function(i) mean(rexp(n))))

   user  system elapsed 
  3.580   0.000   3.665

An alternative is to unlist the result of lapply.

set.seed(123)
str(unlist(lapply(1:N, function(i) mean(rexp(n)))))

 num [1:100000] 0.704 0.584 1.505 1.083 1.489 ...

for which the timing is

system.time(unlist(lapply(1:N, function(i) mean(rexp(n)))))

   user  system elapsed 
  3.430   0.000   3.428

A slight variant on sapply is vapply which can be used if you know you will be creating a vector. It takes a third argument which is a template vector that indicates the desired type of the response. It doesn't need to be as large are the response, it just needs to be the right type of vector (technically, I should say "mode" instead of "type" but you get the idea). Typical values are 1 for a numeric vector, 1L for an integer vector, and the empty string, "", for character strings but more complex structures can be used. See the examples in the help page ?vapply

set.seed(123)
str(vapply(1:N, function(i) mean(rexp(n)), 1))

 num [1:100000] 0.704 0.584 1.505 1.083 1.489 ...

system.time(vapply(1:N, function(i) mean(rexp(n)), 1))

   user  system elapsed 
  3.230   0.000   3.253

Notice that the timings for the sapply and lapply methods are similar to those for replicate. This is not a coincidence. The replicate function

replicate

function (n, expr, simplify = TRUE) 
sapply(integer(n), eval.parent(substitute(function(...) expr)), 
    simplify = simplify)
<environment: namespace:base>

ends up being a call to sapply. There are a few subtleties in there about setting up the counter vector and changing the expression into an anonymous function but basically it comes down to a call to sapply.

We won't bother discussing how the expression is converted to an anonymous function but it is interesting to consider why the number of replications, what we call N but is called n here, is converted to a vector of length n by integer(n), whereas we used 1:N. The integer(n) call creates an integer vector of length n filled with zeros, whereas 1:N creates an integer vector of length N counting from 1 to N — most of the time. Because we are not actually using the values it doesn't matter what the contents of the vector are as long as it has the correct length.

As for that "most of the time" comment, the exception is when n=0. The construction 1:0 produces a vector of length 2 whereas integer(0) produces a vector of length 0. It may seem bizarre to consider what the result should be when you replicate evaluation of an expression 0 times but the defensive programmer is always cautious of the "edge cases". So the construction 1:n is dangerous because it doesn't behave as expected when n=0. If you want to get an integer vector containing the values from 1 up to n with correct behaviour for n=0 use seq_len(n)