Simulating in R- how can I make this faster?

Question

I am simulating something like Jim Berger's applet.

The simulation works like this: I will generate a sample x of size n either from the null distribution N(0,1) or from the alternative distribution N(theta, 1). I will assume that the the prior probability of the null is some proportion prop (so the prior of the alternative is 1-prop) and that the distribution of theta in the alternative is N(0,2) (I could change all this parameters, but this is just to start).

I want to get a large number of pvalues arround a certain range (like 2000 pvalues between 0.049 and 0.05, in the simulation this would be equivalent to z stats arround 1.96 and 1.97) from the simulation scenario described above, and to see how many came from the null and how many came from the alternative.

So far I came up with a solution like this:

berger <- function(prop, n){
  z=0
  while(z<=1.96|z>=1.97){
    u <- runif(1)
    if(u<prop){
      H0 <- TRUE
      x<-rnorm(n, 0, 1)
    }else{
      H0 <- FALSE
      theta <- rnorm(1, 0, 2)
      x <- rnorm(n, theta, 1)
    }
    z <- sqrt(n)*abs(mean(x))
  }
  return(H0)
}

results<-replicate(2000, berger(0.1, 100))
sum(results)/length(results) ## approximately 25%

It takes about 3,5 minutes. Is it possible to speed this up? How? Every answer is welcome, including integration with C.

Update: Parallelization can speed that up a little bit. But, I have tried the same code in Julia, and it takes only 14 seconds without any parallelization (code below).

Update 2: With Rcpp and parallelization it is possible to reduce the simulation to 8 seconds. See the new answer.

function berger(prop, n)
       z = 0 
       h0 = 0
       while z<1.96 || z > 1.97

              u = rand()

              if u < prop
                     h0 = true;
                     x = randn(n)             
              else
                     h0 = false
                     theta = randn()*2
                     x = randn(n) + theta
              end

              z = sqrt(n)*abs(mean(x))
       end

       h0
end

results = [0]

for i in 1:2000
       push!(results, berger(0.1, 100))
end

sum(results)/length(results)

Answer 1

There might be ways to make this function a little faster (through parallelization for example), but you aren't going to get orders of magnitude difference (edit: in R). The key problem is that you are making roughly 400 million draws from the normal distribution.

This is a function that returns the average number of runs through the while your function takes:

f<-function(prop,n){
  i<-0
  z<-0
  while(z<=1.96|z>=1.97){
    i<-i+1
    u <- runif(1)
    if(u<prop){
      H0 <- TRUE
      x<-rnorm(n, 0, 1)
    }else{
      H0 <- FALSE
      theta <- rnorm(1, 0, 2)
      x <- rnorm(n, theta, 1)
    }
    z <- sqrt(n)*abs(mean(x))
  }
  return(i)
}

Now we can calculate how many times your function runs:

set.seed(1)
runs<-replicate(200,f(prop=0.1, n=100))
mean(runs) # 2034
sd(runs) # 2121

So, to calculate the number of draws from the normal distribution:

# number of replicates
# times normal distributions per replicate
# draws from each distribution
2000*mean(runs)*100
# 406,853,000 normal distribution draws

The rnorm function calls a compiled C function, and is likely to be near optimal speed. You can test the "lower bound" of making this many draws on your own machine:

system.time(rnorm(406853000))
# My machine:
#   user  system elapsed 
#  53.78    2.39   56.62 

By comparison your function runs roughly four times slower:

system.time(replicate(2000,berger(prop=0.1,n=100)))
#    user  system elapsed 
#  210.40    0.03  211.12 

So, your function really isn't that slow when you think about it, especially when you consider that there is overhead on each call to rnorm. If it is very critical that you improve the speed of this function, and you have a few cores, you can easily parallelize it in R:

library(parallel)
mclapply(1:2000,function(x) berger(prop=0.1,n=100))

Other than that, you could write a super-optimized function in C and save a few minutes, but it might not be worth it.

Answer 2

It is actually simple to speed this up using Rcpp. Combining Rcpp with parellelization, I was able to reduce the time to 8 seconds.

The .cpp file is something like this (using Rcpp "sugars" turns this task pretty easy - since this was the first time I used Rcpp, maybe this code is not optimal, but it did the job!):

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]


int RcppBerger(double prop, int n) {

  double z=0,theta=0, u=0;
  int h = 0;
  NumericVector x;
    while (z<1.96 || z > 1.97){
      u = R::runif(0, 1);
      if(u < prop){
        h = 1;
        x = rnorm(n);
        }else{
          h = 0;
          theta = R::rnorm(0, 2);
          x = rnorm(n, theta, 1);
          }
          z = sqrt(n)*mean(x);
          if(z<0){z = -1*z;};
    }
  return h;
}

Then, without parallelization, you can just use the sourceCpp function and the RcppBerger will be available in the workspace:

library(Rcpp)
sourceCpp("RcppBerger.cpp")
results<-replicate(2000, RcppBerger(0.1, 100))
sum(results)/length(results) ## approximately 25%

This already reduces the time from 3,5 minutes to arround 40 seconds. After that we can parallelize.

In Windows, this is a little tricky, it seems that you have to create a package first. But Rcpp provides a nice function to do that Rcpp.package.skeleton. Just put the source file in it and it will create all the necessary documents and folders:

Rcpp.package.skeleton("RcppBerger", cpp_files = "RcppBerger.cpp")

Then, after installing the package, you can parallelize with foreach and doParallel:

library(foreach)
library(doParallel)
library(RcppBerger)
registerDoParallel(cores=8)
results<- foreach(1:2000, .packages="RcppBerger") %dopar% RcppBerger(0.1, 100)

Now the simulation takes only 8 seconds.