rcpp

R 中的常用命令

rm(list = ls()) #清空变量

xVector1 = c(1,2,3,4,5) # num格式的向量
xVector1

xVector2 = 1:7 # int格式的向量
xVector2


yMatrix = matrix(c(1,2,3,4,5,6), 
                 nrow = 2,
                 ncol = 3
                 )
yMatrix # R中的矩阵按列生成


aList = list(aVector = c(1,3,5),
             aMatrix = matrix(1:4,nrow = 2))
aList
aList$aVector # $ 符号来访问list容器中的元素
aList$aMatrix

library(microbenchmark)

mean_forloop = function(x){
    n = length(x)
    s = 0
    for (i in 1:n) {
        s = s + x[i]
    }
    return(s/n)
}

mean_direct = function(x){
    return(mean(x))
}

p = 1e4
xData = runif(p)
testResult = microbenchmark(mean_forloop(xData),
                            mean_direct(xData),
                            times = 50)
testResult

例1

使用 RcppArmadillo 在 R 中写 c++ 的项目

#include <RcppArmadillo.h> 
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;

// [[Rcpp::export]]
vec timesTwo(vec x){
    vec y;
    y = 2*x;
    return y;
}

• #include <RcppArmadillo.h> 使用 RcppArmadillo 所需的头文件
• // [[Rcpp::depends(RcppArmadillo)]] 提供了链接的信息
• // [[Rcpp::export]] 表明紧随其后的函数将输出到 R 当中，可以被 R 调用

例2

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;

// [[Rcpp::export]]
mat matrixTimesTwo(mat x){
    mat y;
    y = 2*x;
    return y;
}

例3 平方

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
#include <cmath>
using namespace arma;

// [[Rcpp::export]]
vec vectorSquared (vec x){
    int i, vSize = x.n_elem;
    vec xSq(vSize);
    for(i=0; i<vSize; i++){
        xSq(i) = pow(x(i),2);
    }
    return xSq;
}

// [[Rcpp::export]]
mat matrixSquared (mat x){
    int i,j,m,n;
    m = x.n_rows;
    n = x.n_cols;
    mat xSq(m, n);
    for(i=0; i<m; i++)
        for(j=0; j<n; j++)
            xSq(i,j) = pow(x(i,j), 2);
    return xSq;
}

例4 卷积速度

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;

// [[Rcpp::export]]
double convolve_c(vec A, vec B){
    size_t vSize, i;
    double result = 0;
    
    vSize = A.size();
    for(i=0; i<vSize; i++){
        result += A(i) * B(vSize - i - 1);
    }
    return result;
}

convolve_r = function(A, B){
    result = convolve(A,B)
    return(result)
}

p = 1e4
x <- runif(p)
y <- runif(p)

testConv = microbenchmark(convolve_c(x,y),
                          convolve_r(x,y),
                          # simple_r(x,y),
                          times = 50)
testConv

例5 局部平滑

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;

int mirrorIndex(int fetchI, int length){
  if(fetchI < 0)
    fetchI = -fetchI - 1;
  if(fetchI >= length){
    fetchI = length - (fetchI - length) - 1;
  }
  return fetchI;
}

// [[Rcpp::export]]
vec localSmoothing(vec y, vec weights){
  int i, j, ySize = y.n_elem;
  vec smoothedY(ySize);
  double s;
  int fetchI;
  for(i = 0; i < ySize; i++){
    s = 0;
    for(j = -4; j < 5; j++){
      fetchI = i + j;
      fetchI = mirrorIndex(fetchI, ySize);
      s += weights(j+4) * y(fetchI);
    }
    smoothedY(i) = s;
  }
  
  return smoothedY;
}

rm(list = ls()) #清空变量

set.seed(100)
n = 400
x = 2*pi*(1:n)/n
y = sin(x) + rnorm(n, sd = 0.1)

weights = c(1,2,3,4,5,4,3,2,1) /25
fHat = localSmoothing(y, weights)

par(mar = c(4,4,1,1))
plot(x,y, pch = 20)
lines(x, sin(x), lwd = 2, lty = 2, col = "blue")
lines(x,fHat, lwd = 2, col = "red")