Rcpp-training.Rmd

---
title: "Rcpp training"
author: Jeanne Clément and Ghislain Vieilledent
date: "`r format(Sys.time(),'%B %e, %Y')`"
output:
  bookdown::html_document2:
    #base_format: rmarkdown::html_vignette
    highlight: tango
    number_sections: true
    toc: true
    toc_float: true
    fig_caption: yes
link-citations: yes
urlcolor: Maroon
---

```{r options, include=FALSE}
library(knitr)
library(kableExtra)
library(bookdown)
library(dplyr)

opts_chunk$set(echo=TRUE, cache=FALSE,
               #results="hide", 
               warning=FALSE,
               message=FALSE, highlight=TRUE,
               fig.show="hide", size="small",
               fig.align="center",
               tidy=FALSE)
options(knitr.kable.NA="-")
```

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/. 

# 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

# RcppGSL example

Function to generate a random sample from gaussian distribution: 

## C++ and R code

```{Rcpp RcppGSL-my_rnorm}
#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 + sigma * gsl_ran_ugaussian(s); // 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)
*/
```

## Representation of results 

```{r r-my_rnorm, fig.show="asis", out.width="70%"}
hist(beta)
```

# RcppArmadillo example

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

## C++ code

```{Rcpp RcppArma-example, echo=T}
#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;
}

```

## R code 

```{r r-RcppArma-example}
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)
var(X_scaled)
```

# Distance computation 

## 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}.$$

## R function

```{r R_distmat}

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)
}
```

## C++ function 

```{Rcpp arma_distmat}
#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;
}
```

## Data simulation 

**Simulation of 500 points dispersed all over the space**

```{r points-representation, fig.show="asis", out.width="60%", out.height="60%", fig.align="left"}
# 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**

```{r coordinates, echo=FALSE}
knitr::kable(coords[1:5,],col.names = c(" x "," y "), digits=1, booktabs=TRUE) %>%
		kableExtra::kable_styling(latex_options=c("HOLD_position","striped"),
		                          full_width=FALSE,position = "left")
```

## Comparison of compilation times

```{r dist-benchmark}
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"))
```

```{r dist-benchmark-results, echo=FALSE}
library(dplyr)
knitr::kable(Benchmark,digits=3, booktabs=TRUE) %>%
		kableExtra::kable_styling(latex_options=c("HOLD_position","striped"),
		                          full_width=FALSE,position = "left")
```

## Results 

**Distance matrix for the 5 points**

```{r distance-matrix, echo=F}
distmat <- arma_distmat(coords)
knitr::kable(distmat[1:5,1:5], digits=1, booktabs=TRUE) %>%
		kableExtra::kable_styling(latex_options=c("HOLD_position","striped"),
		                          full_width=FALSE,position = "left")
```

# Simple linear regression 

## 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}.$$

## R function

```{r r-fastLm}
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))
}
```


## C++ function

```{Rcpp fastLm}
#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);
}
```

## 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. 

```{r lin-trees}
library(datasets)

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

## Comparison of compilation times

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

```{r lin-benchmark-results, echo=F}
knitr::kable(Benchmark,digits=4, booktabs=TRUE) %>%
		kableExtra::kable_styling(latex_options=c("HOLD_position","striped"),
		                          full_width=FALSE,position = "left")
```

## Representation of results 

```{r lin-results,fig.show="asis", out.width="50%",out.height="50%"}

# 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")
```

# Log-likelihood computation 

## 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}$$

## 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. 

## R function

```{r R_logL}
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)
}
```

## C++ function

Implement the same algorithm in C++ using RcppArmadillo.