In hudie-lab/StatComp: Example functions for 'Statistical Computing' course

Overview

StatComp is a simple R package developed to compare the performance of R and R++ (implemented through the R package Rcpp) for the 'Statistical Computing' course. Two functions are considered, namely, gibbs (generate random nubers using Gibbs sampler) and vacc (predicting the response using three variables: age, gender, and ily). For each function, both R and Rcpp versions are produced. Namely gibbsR and vaccR for R and 'gibbC' and vaccC for C++.

The R package 'microbenchmark' can be used to benchmark the above R and C++ functions.

Benchmarking vaccR and vaccC

The source R code for vaccR is as follows:

function (age, female, ily) 
{
    p <- 0.25 + 0.3 * 1/(1 - exp(0.04 * age)) + 0.1 * ily
    p <- p * ifelse(female, 1.25, 0.75)
    p <- pmax(0, p)
    p <- pmin(1, p)
    p
}

The above code involves 'ifelse', 'pmax', and 'pmin', which are known to be slow. On the other hand, the following Rcpp code is much faster.

double vacc3a(double age, bool female, bool ily){
  double p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily;
  p = p * (female ? 1.25 : 0.75);
  p = std::max(p, 0.0);
  p = std::min(p, 1.0);
  return p;
}
NumericVector vaccC(NumericVector age, LogicalVector female,
                    LogicalVector ily) {
  int n = age.size();
  NumericVector out(n);
  for(int i = 0; i < n; ++i) {
    out[i] = vacc3a(age[i], female[i], ily[i]);
  }
  return out;
}

In order to empirically benchmark vaccR and vaccC, one generates 1,000 repicates of (age, female, ily), and save it in data{StatComp}. To load the file, one simply calls data(data). The R code for benchmark vaccR and vaccC is as follows.

library(StatComp)
library(microbenchmark)
data(data)
attach(data)
tm2 <- microbenchmark(
  vR = vaccR(age,female,ily),
  vC = vaccC(age,female,ily)
)
knitr::kable(summary(tm2)[,c(1,3,5,6)])

The above results show an evident computational speed gain of C++ against R.

Benchmarking gibbsR and gibbsC

The source R code for vaccR is as follows:

gibbsR <- function(N, thin) {
  mat <- matrix(nrow = N, ncol = 2)
  x <- y <- 0
  for (i in 1:N) {
    for (j in 1:thin) {
      x <- rgamma(1, 3, y * y + 4)
      y <- rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))
    }
    mat[i, ] <- c(x, y)
  }
  mat
}

The above R code involves two loops, which could be very slow even for R-3.01 or any higher version. The corresponding C++ code is as follows.

NumericMatrix gibbsC(int N, int thin) {
  NumericMatrix mat(N, 2);
  double x = 0, y = 0;
  for(int i = 0; i < N; i++) {
    for(int j = 0; j < thin; j++) {
      x = rgamma(1, 3, 1 / (y * y + 4))[0];
      y = rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))[0];
    }
    mat(i, 0) = x;
    mat(i, 1) = y;
  }
  return(mat);
}

The R code for benchmarking gibbsR and gibbsC is as follows.

tm2 <- microbenchmark(
  vR = gibbsR(1e4, 10),
  vC = gibbsC(1e4, 10)
)
knitr::kable(summary(tm2)[,c(1,3,5,6)])

The results again show an evident computational speed gain of C++ against R.

hudie-lab/StatComp documentation built on Jan. 3, 2020, 12:15 a.m.