In hudie-lab/SC19036: The R Package For StatComp

HW THEME

• Rcpp

QUESTION

• Rewrite a Rcpp function to implement a random walk Metropolis sampler for generating the standard Laplace distribution.
• Compare the results by the two functions: Rcpp function and R function.
  • Compare the generated random numbers by the two functions using qqplot.
  • Campare the computation time of the two functions with microbenchmark.

ANSWER

function of R

dl <- function (x) 0.5*exp(-abs(x))
rw.Metropolis <- function( sigma, x0, N) { 
  x <- numeric(N)
  x[1] <- x0
  u <- runif(N)
  k <- 0
  for (i in 2:N) {
    y <- rnorm(1, x[i-1], sigma)
    if (u[i]*dl(x[i-1]) <= dl(y)) {x[i]<-y} 
    else{
    x[i] <- x[i-1]
    k <- k + 1 }
  }
  return(list(x=x, k=k)) 
}

function of C++

library(Rcpp)
cppFunction('NumericVector laprn (double sigma,double x0,int N) {
    std::unordered_set<int> seen;
    NumericVector x(N);
    NumericVector y(N);
    x[0] = x0;
    NumericVector u = runif(N);
    int k = 0;
    for (int i =1; i < N;i++) {
        NumericVector y = rnorm(1, x[i-1], sigma);
        if (u[i-1] <= (0.5*exp(-abs(y[0]))) / (0.5*exp(-abs(x[i-1])))){
            x[i] = y[0];
        }
        else{
          x[i] = x[i-1];
          k = k + 1 ;
          }
    }
    return(x);
}')

comparison of the computation time

    # Can create source file in Rstudio
    library(microbenchmark)
    x0 <- 25
    N <- 2000
    sigma <- 0.05
    (ts <- microbenchmark(rnR=rw.Metropolis(sigma,x0,N),
                         rnC=laprn(sigma,x0,N)))

comparison by qqplot

rnR=rw.Metropolis(sigma,x0,N)[[1]]
rnC=laprn(sigma,x0,N)
qqplot(rnR,rnC)

Footnote : the X-axis stands for the random number by R, the Y-axis by C++. And here we generate r N numbers for each method and set sigma as r sigma.

Comments : the QQplot shows that the random numbers by Rcpp run quicker to the stable chains and the method used by R runs slower and is more difficult to arrive the stable situation.

hudie-lab/SC19036 documentation built on Jan. 2, 2020, 8:45 p.m.