RcppArmadillo(Rcpp) を使用してコードを高速に実行する方法を学んでいます。実際には、特定のパッケージからいくつかの関数を呼び出したい場合がよくあります。以下はほんの一例です。lasso (または scad、mcp) のしきい値を計算したい。R では、thresh_est
関数 inを使用できlibrary(HSDiC)
ます。
library(HSDiC) # acquire the thresh_est() function
# usage: thresh_est(z, lambda, tau, a = 3, penalty = c("MCP", "SCAD", "lasso"))
# help: ?thresh_est
z = seq(-5,5,length=500)
thresh = thresh_est(z,lambda=1,tau=1,penalty="lasso")
# thresh = thresh_est(z,lambda=1,tau=1,a=3,penalty="MCP")
# thresh = thresh_est(z,lambda=1,tau=1,a=3.7,penalty="SCAD")
plot(z,thresh)
次に、RcppArmadillo(Rcpp) を介して上記を実現しようとします。teuderと coatless およびcoatlessの回答によると 、コード ( thresh_arma.cppとして保存されます) を以下に記述します。
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp; //with this, Rcpp::Function can be Function for short, etc.
using namespace arma; // with this, arma::vec can be vec for short, etc.
// using namespace std; //std::string
// [[Rcpp::export]]
vec thresh_arma(vec z, double lambda, double a=3, double tau=1, std::string penalty) {
// Obtaining namespace of the package
Environment pkg = namespace_env("HSDiC");
// Environment HSDiC("package:HSDiC");
// Picking up the function from the package
Function f = pkg["thresh_est"];
// Function f = HSDiC["thresh_est"];
vec thresh;
// In R: thresh_est(z, lambda, tau, a = 3, penalty = c("MCP", "SCAD", "lasso"))
if (penalty == "lasso") {
thresh = f(_["z"]=z,_["lambda"]=lambda,_["a"]=a,_["tau"]=tau,_["penalty"]="lasso");
} else if (penalty == "scad") {
thresh = f(_["z"]=z,_["lambda"]=lambda,_["a"]=a,_["tau"]=tau,_["penalty"]="SCAD");
} else if (penalty == "mcp") {
thresh = f(_["z"]=z,_["lambda"]=lambda,_["a"]=a,_["tau"]=tau,_["penalty"]="MCP");
}
return thresh;
}
次に、コンパイルを次のように行います
library(Rcpp)
library(RcppArmadillo)
sourceCpp("thresh_arma.cpp")
# z = seq(-5,5,length=500)
# thresh = thresh_arma(z,lambda=1,tau=1,penalty="lasso")
# # thresh = thresh_arma(z,lambda=1,a=3,tau=1,penalty="mcp")
# # thresh = thresh_arma(z,lambda=1,a=3.7,tau=1,penalty="scad")
# plot(z,thresh)
ただし、コンパイルは失敗し、理由はわかりません。