Rcpp training

A short presentation of Rcpp* packages is available at: https://jeanneclement.github.io/Rcpp-training/pres_Rcpp.pdf.

This tutorial is available at the following webpage: https://jeanneclement.github.io/Rcpp-training/.

All the source code is available on GitHub: https://github.com/JeanneClement/Rcpp-training/.

1 Packages to install

• Rcpp
• RcppGSL require GSL installed (~$ sudo apt-get install libgsl23 libgsl-dev with Debian/Ubuntu)
• RcppArmadillo
• datasets
• rdist
• rbenchmark
• rmarkdown
• knitr
• kableExtra
• magrittr

2 RcppGSL example

Function to generate a random sample from gaussian distribution:

2.1 C++ and R code

#include <Rcpp.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

// [[Rcpp::depends(RcppGSL)]]

// [[Rcpp::export]]
Rcpp::NumericVector my_rnorm(const int& nsamp, const double& mu,
                             const double& sigma, const int& seed) {
  
  // Initialize random seed for random number generator 
  gsl_rng *s = gsl_rng_alloc(gsl_rng_mt19937);
  gsl_rng_set(s, seed);
  
  // Initialize the vector to store results
  Rcpp::NumericVector beta(nsamp);
  
  // Generate a sample of nsamp values from gaussian distribution 
  for (int i = 0; i < nsamp; i++) {
    beta[i] = mu + gsl_ran_gaussian(s, sigma); // Random draw
  }
  
  // Free memory
  gsl_rng_free(s);
  
  return beta;
}

// Possibility to add R code in .cpp file or Rcpp chunk
// Rcpp::sourceCpp() compile the C++ function then execute R code 

/*** R 

library(Rcpp)
library(RcppGSL)

## Draw random sample 
beta <- my_rnorm(nsamp=100, mu=5, sigma=2, seed=123)

*/

2.2 Representation of results

hist(beta)

3 RcppArmadillo example

Function to center and scale the columns of a numeric matrix:

3.1 C++ code

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat arma_scale(const arma::mat& X) {
  
  // Define dimensions of the matrix 
  int n = X.n_rows, k = X.n_cols;
  
  // Initialize the matrix to store results
  arma::mat X_scaled(n,k);
  
  // means 
  arma::rowvec col_means = arma::mean(X,0); 
  // standard deviations
  arma::rowvec col_sd  = arma::stddev(X,0); 
  
  // Scale and center each column of the matrix 
  for (int  p= 0; p < k; p++) {
    X_scaled.col(p) = (X.col(p)-col_means(p))/col_sd(p);
  }
  
  return X_scaled;
}

3.2 R code

library(Rcpp)
library(RcppArmadillo)

# Center and reduce a matrix 
X <- matrix(rnorm(50),ncol=5)
X_scaled <- arma_scale(X)
## Check if the resulting matrix is scaled and centered
colMeans(X_scaled)

## [1] -3.469447e-17 -1.137979e-16  2.268151e-17  5.551115e-18  2.428613e-17

var(X_scaled)

##             [,1]         [,2]        [,3]       [,4]         [,5]
## [1,]  1.00000000  0.632458512  0.05379625 -0.5393759  0.163506516
## [2,]  0.63245851  1.000000000 -0.05137871 -0.4164714  0.006563847
## [3,]  0.05379625 -0.051378708  1.00000000 -0.4441156  0.762666664
## [4,] -0.53937594 -0.416471383 -0.44411557  1.0000000 -0.480016645
## [5,]  0.16350652  0.006563847  0.76266666 -0.4800166  1.000000000

4 Distance computation

4.1 Mathematical definition

The Euclidean distance between two points whose coordinates are \(A=(x_A, y_A)\) and \(B=(x_B,y_B)\) is given by \[\sqrt{(x_B-x_A)^2 + (y_B-y_A)^2}.\]

4.2 R function

R_distmat <- function(X) {
  
  # Number of points
  np <- nrow(X)
  
  # Initialize with zeros the matrix to store results
  distmat <- matrix(0,nrow=np,ncol=np)
  
  # Loop on points
  for (i in 1:(np-1)) {
    p0 <- X[i,] # fix a point 
    
    # Loop to calculate the distances between this point and the next ones  
    for (j in (i+1):np){
    p1 <- X[j,]
    diff <- p0-p1 # (x0-x1,y0-y1)
    squared_diff <- t(diff)%*% diff # (x0-x1)² + (y0-y1)²
    
    # Fill the distance matrix with the square root of precedent value
    distmat[j,i] <- distmat[i,j] <- sqrt(squared_diff)
    }
  }
 return(distmat)
}

4.3 C++ function

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
arma::mat arma_distmat(const arma::mat& X) { 

  // Number of points
  int np = X.n_rows;
  
  // Initialize with zeros the matrix to store results
  arma::mat distmat; distmat.zeros(np, np);
  
  // Loop on all points 
  for (int i = 0; i < np; i++) {
    arma::vec p0 = X.row(i).t(); // fix a point 
    
    // Loop to calculate the distances between this point and the next ones  
    for (int j = i + 1; j < np; j++) {
      arma::vec p1 = X.row(j).t();
      arma::vec diff = p0 - p1; // (x0-x1,y0-y1)
      double squared_diff = as_scalar(diff.t() * diff); // (x0-x1)² + (y0-y1)²
      // Fill the distance matrix with the square root of precedent value
      distmat(j, i) = distmat(i, j) = sqrt(squared_diff);
    }
  }
  return distmat;
}

4.4 Data simulation

Simulation of 500 points dispersed all over the space

# Data simulation 
coords <- matrix(runif(500,0,100),ncol=2)

# Spatial representation of points  
plot(coords[,1],coords[,2],pch=4,xlab="x", ylab="y",
     main="Spatial repartition of points", col="forestgreen")

Coordinates of 5 points

x	y
58.5	55.5
15.2	22.0
5.8	11.6
65.9	83.0
88.9	95.1

4.5 Comparison of compilation times

library(rdist)

# Benchmark
library(rbenchmark)
Benchmark <- benchmark(
  "arma_distmat" = {arma_distmat(coords)},
  "cdist" = {cdist(coords,coords)},
  "R_distmat" = {R_distmat(coords)},
  replications=10,
  columns = c("test", "elapsed", "relative"))

test	elapsed	relative
arma_distmat	0.009	1.000
cdist	0.047	5.222
R_distmat	1.116	124.000

4.6 Results

Distance matrix for the 5 points

0.0	54.7	68.5	28.5	49.9
54.7	0.0	14.0	79.3	103.8
68.5	14.0	0.0	93.3	117.7
28.5	79.3	93.3	0.0	26.0
49.9	103.8	117.7	26.0	0.0

5 Simple linear regression

5.1 Mathematical definition

We have \(n\) observations of a response variable \(Y=(y_i)_{i=1,\ldots,n}\) and \(p\) explanatory variables \((X_1,\ldots,X_p)\),
such as \(X_1=(x_{11},\ldots,x_{i1},\ldots,x_{n1})'\).

We want to estimate coefficients of the linear regression \(\beta=(\beta_0,\beta_1,\ldots,\beta_p)'\) such as : \[ y_i = \beta_0 + \beta_1x_{i1}+\ldots+\beta_px_{ip}+\epsilon_i,\] where \(\epsilon_i \sim \mathcal{N}(0,\sigma^2) \ iid\).

Then \(y_i \sim \mathcal{N}(\beta_0 + \beta_1x_{i1}+\ldots+\beta_px_{ip}, \ \sigma^2)\) \(iid\) for \(i=1,\ldots,n\).

This gives in matrix writing :

\[Y = X\beta + \epsilon \] where \(X=(\mathbb{1}_n,X_1,\ldots,X_p)\) and \(\epsilon=(\epsilon_i)_{i=1,\ldots,n}\).

According to the Ordinary Least Squares (OLS) method, \(\beta\) is estimated by : \[\widehat{\beta}=\left(X'X\right)^{-1}X'y.\]

We define the residuals \[\widehat{\epsilon} = Y - X\widehat{\beta}=Y-\widehat{Y}.\]

Then residual variance is given by \[\widehat{\sigma^2}=\dfrac{\sum\limits_{i=1}^n \widehat{\epsilon_i}^2}{n-p-1}.\]

Finally the variance-covariance matrix of coefficients is estimated by \[V_{\widehat{\beta}}=\widehat{\sigma^2}\left(X'X\right)^{-1}.\]

5.2 R function

R_fastLm <- function(X, y) {
n <- nrow(X)
p <- ncol(X)

# fit model y ~ X       
coef <- solve(t(X) %*% X) %*% t(X) %*% y

# residuals
res <- y - X %*% coef

# std.errors of coefficients
s2 <- sum(t(res) %*% res)/(n - p);
std_err <- sqrt(s2 %*% diag(solve(t(X) %*% X)))

return(list("coefficients" = coef,
            "stderr" = std_err,
            "residuals" = res,
            "sigma2" = s2,
            "df.residual" = n - p))
}

5.3 C++ function

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
Rcpp::List arma_fastLm(const arma::mat& X, const arma::colvec& y) {
    int n = X.n_rows;
    int p = X.n_cols;
        
    arma::colvec coef = arma::inv(X.t()*X)*X.t()*y; // fit model y ~ X
    arma::colvec res  = y - X*coef; // residuals
    
    // std.errors of coefficients
    double s2 = arma::sum(res.t()*res)/(n-p);
    arma::colvec std_err = arma::sqrt(s2*arma::diagvec(arma::inv(X.t()*X)));
    
    return Rcpp::List::create(Rcpp::Named("coefficients") = coef,
                        Rcpp::Named("stderr") = std_err,
                        Rcpp::Named("residuals")  = res,
                        Rcpp::Named("sigma2")  = s2,
                        Rcpp::Named("df.residual")  = n - p);
}

5.4 Linear regression on trees data set

This data set provides measurements of the girth, height and volume of timber in 31 felled black cherry trees.

Volume will be considered as the variable to be explained and girths and heights as the explanatory variables. To have linear relationship between this variables, we will consider their log.

library(datasets)

# Trees data-set
y <- log(trees$Volume)
X <- cbind(1, log(trees$Girth), log(trees$Height))

5.5 Comparison of compilation times

# Benchmark
library(rbenchmark)
Benchmark <- benchmark(
  "arma_fastLm" = {arma_fastLm(X,y)},
  "R_fastLm" = {R_fastLm(X,y)},
  replications=100,
  columns = c("test", "elapsed", "relative"))

test	elapsed	relative
arma_fastLm	0.002	1
R_fastLm	0.008	4

5.6 Representation of results

# Trees data-set
y <- log(trees$Volume)
X <- cbind(1, log(trees$Girth), log(trees$Height))

# fit model y ~ X
arma_mod <- arma_fastLm(X,y)
R_mod <- R_fastLm(X,y)

# Comparison of residuals obtained with R and C++ functions
plot(arma_mod$residuals, R_mod$residuals,
     xlab="residuals with C++", ylab="residuals with R",
     pch=4, main="Residuals")
abline(a=0,b=1, col="red")

# Representation of the linear regression
plot(X%*%arma_mod$coefficients, y, xlab="Xbeta",
     ylab="y", pch=4, main="Linear regression")
abline(a=0,b=1, col="red")

6 Log-likelihood computation

6.1 Mathematical definition

The likelihood function expresses the plausibilities of different parameter values for a given sample of data. The maximum of this function, if it exists, correspond to the combination of model parameter values that maximize the probability of drawing the sample actually obtained.

The likelihood corresponding to the previous simple linear model is given by : \[\begin{aligned} L(\beta,\sigma^2) &= \prod \limits_{i=1}^n p(y_i \ | \ \beta, \sigma^2) \\ &=\prod \limits_{i=1}^n \frac{1}{\sigma \sqrt{2\pi}}\exp{\left(-\frac{1}{2\sigma^2}(y_i-\beta_0-x_{i1}\beta_1-\ldots-x_{ip}\beta_p)^2 \right)} \end{aligned}\]

Then the log-likelihood is : \[\begin{aligned} l(\beta,\sigma^2) &= \log{\left(L(\beta,\sigma^2)\right)} \\ &= \sum \limits_{i=1}^n \ \log \left(\frac{1}{\sigma \sqrt{2\pi}}\right) -\frac{1}{2\sigma^2}(y_i-\beta_0-x_{i1}\beta_1-\ldots-x_{ip}\beta_p)^2\\ &= -n\log(\sigma)-n\frac{\log(2\pi)}{2} -\frac{1}{2\sigma^2} \sum \limits_{i=1}^n (y_i-\beta_0-x_{i1}\beta_1-\ldots-x_{ip}\beta_p)^2 \end{aligned}\]

6.2 Exercise

1. Implement a function in C++ using RcppArmadillo to compute the log-likelihood corresponding to the previous simple linear model taking as argument the explicative variables \(X\), the response variable \(Y\), the estimated regression coefficients \(\widehat{\beta}\) and the computed residual variance \(\widehat{\sigma^2}\).

2. Use trees data set to fit three distincts linear models with the log of the Volume as response variable and different choices of explanatory variables.
   

3. Compute the log-likelihood of each model.

4. Deduce from the values obtained which model best reflects the data.

6.3 R function

R_logL <- function(X,y,beta,sigma2) {
  
  n <- length(y)
  p <- ncol(X)
  
  logL <- 0
  for (i in 1:n) {
    Xbeta_part <- 0
    for (j in 1:p) {
      Xbeta_part <- Xbeta_part + X[i,j]*beta[j]
    }
    logL <- logL + dnorm(x=y[i], mean=Xbeta_part, sd=sqrt(sigma2), log=T)
  }
  
  return(logL)
}

6.4 C++ function

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
double arma_logL(const arma::mat& X, const arma::vec& y,
                 const arma::vec& beta, const double& sigma2) {
  
  int n = y.n_elem;
  int p = X.n_cols;
  
  double logL = 0;
  for (int i = 0; i < n; i++) {
      double Xbeta_part = 0;
    for (int j = 0; j < p; j++) {
      Xbeta_part += X(i,j)*beta(j);
    }
    logL += R::dnorm(y[i], Xbeta_part, std::sqrt(sigma2), 1);
  }
  
  return logL;
}

6.5 Linear regressions on trees data set

# Trees data-set
y <- log(trees$Volume)

# fit models y ~ X
## model 1 : log(Volume) ~ log(Girth)
X1 <- cbind(1, log(trees$Girth))
arma_mod1 <- arma_fastLm(X1,y)
R_mod1 <- R_fastLm(X1,y)

## model 2 : log(Volume) ~ log(Height)
X2 <- cbind(1, log(trees$Height))
arma_mod2 <- arma_fastLm(X2,y)
R_mod2 <- R_fastLm(X2,y)

## model 3 : log(Volume) ~ log(Height) + log(Girth)
X3 <- cbind(1, log(trees$Girth), log(trees$Height))
arma_mod3 <- arma_fastLm(X3,y)
R_mod3 <- R_fastLm(X3,y)

6.6 Comparison of compilation times

We compare the speed of C++ and R functions used to compute de log-likelihood corresponding to model 3.

# Benchmark
Benchmark <- benchmark(
  "arma_logL" = {arma_logL(X3,y,arma_mod3$coefficients,arma_mod3$sigma2)},
  "R_logL" = {R_logL(X3,y,R_mod3$coefficients,R_mod3$sigma2)},
  replications=100,
  columns = c("test", "elapsed", "relative"))

test	elapsed	relative
arma_logL	0.009	1.8
R_logL	0.005	1.0

6.7 Model performance based on likelihood

# log-likelihood computation
## model 1 : log(Volume) ~ log(Girth)
R_logL1 <- R_logL(X1,y,R_mod1$coefficients, R_mod1$sigma2)
arma_logL1 <- arma_logL(X1,y,arma_mod1$coefficients,arma_mod1$sigma2)

## model 2 : log(Volume) ~ log(Height) 
arma_logL2 <- arma_logL(X2,y,arma_mod2$coefficients,arma_mod2$sigma2)
R_logL2 <- R_logL(X2,y,R_mod2$coefficients,R_mod2$sigma2)

## model 3 : log(Volume) ~ log(Height) + log(Girth)
R_logL3 <- R_logL(X3,y,R_mod3$coefficients,R_mod3$sigma2)
arma_logL3 <- arma_logL(X3,y,arma_mod3$coefficients,arma_mod3$sigma2)

# Results 
logL <- data.frame(arma_logL=c(arma_logL1,arma_logL2,arma_logL3),
                   R_logL=c(R_logL1,R_logL2,R_logL3))
rownames(logL) <- c("model 1","model 2","model 3")
knitr::kable(logL,digits=3, booktabs=TRUE, row.names = T) %>%
        kableExtra::kable_styling(latex_options=c("HOLD_position","striped"),
                                  full_width=FALSE,position = "left")

	arma_logL	R_logL
model 1	24.072	24.072
model 2	-15.149	-15.149
model 3	35.278	35.278

Model 3 has the highest likelihood, it is the the model that best matches the data.