Вызов glmnet от Rcpp (Armadillo) - PullRequest
Новая
       28

Вызов glmnet от Rcpp (Armadillo)

1 голос
/ 10 июня 2019

Я хотел бы извлечь оценки коэффициентов glmnet (после перекрестной проверки) в Rcpp Armadillo, чтобы использовать их в другой функции Armadillo.Я искал похожий вопрос, но не смог найти решение.

Я прилагаю попытку.(Не работает) Даже если бы я получил результат списка cv.glmnet, я бы не смог использовать функцию coef для получения коэффициентов.

R код 

library(glmnet)

set.seed(1)
X = matrix(rnorm(1e3 * 201), 1e3, 201)
beta = -100:100
y = X%*%beta + rnorm(1e3)
cvfit = cv.glmnet(X, y, alpha = 1)
coefs = coef(cvfit, s = "lambda.min")
coefs                                   # get these coefficients from Rcpp

Аргументы cv.glmnet

args(cv.glmnet)
> function (x, y, weights, offset = NULL, lambda = NULL, type.measure = c("mse", "deviance", "class", "auc", "mae"), nfolds = 10, foldid, 
    alignment = c("lambda", "fraction"), grouped = TRUE, keep = FALSE, 
    parallel = FALSE, ...) 
NULL

C ++ код 

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

// [[Rcpp::export]]
Rcpp::List f_cpp(const arma::mat &x, const arma::vec &y, 
                 const arma::vec &weights, 
                 const arma::vec &lambda, double alpha, 
                 int nfolds = 10){

  Rcpp::Environment pkg = Rcpp::Environment::namespace_env("glmnet");

  Rcpp::Function f_R = pkg["cv.glmnet"];

  Rcpp::Nullable<arma::vec> offset = pkg["offset"];
  Rcpp::CharacterVector type_measure = pkg["type.measure"];
  arma::vec foldid = pkg["foldid"];
  Rcpp::CharacterVector alignment = pkg["alignment"];
  bool grouped = pkg["grouped"];
  bool keep = pkg["keep"];
  bool parallel = pkg["parallel"];

  return f_R(x, y, weights, offset, lambda, 
             type_measure, nfolds, foldid, 
             alignment, grouped, keep, parallel, alpha = alpha);
  // coef(f_R(...)) ???
}

1 Ответ

2 голосов
/ 12 июня 2019

Вызов таких функций, как cv.glmnet из C ++, сложен (возможно, даже невозможен), поскольку он использует довольно много возможностей, предлагаемых R, которые делают сигнатуру функции очень гибкой. Однако в R можно определить функцию-обертку, которая использует фактически использованную сигнатуру. Вместо того, чтобы получать эту функцию из (глобальной) среды, я предпочитаю передавать ее как параметр функции: 

library(glmnet)
#> Loading required package: Matrix
#> Loading required package: foreach
#> Loaded glmnet 2.0-16

set.seed(1)
X = matrix(rnorm(1e3 * 201), 1e3, 201)
beta = -100:100
y = X%*%beta + rnorm(1e3)


# set seed since cv.glmnet uses random numbers
set.seed(1)
cvfit = cv.glmnet(X, y, alpha = 1)
coefs = coef(cvfit, s = "lambda.min")

# set seed since cv.glmnet uses random numbers
set.seed(1)
my.glmnet <- function(x, y, alpha) {
    cvfit <- cv.glmnet(x, y, alpha = alpha)
    coef(cvfit, s = "lambda.min")
}
Rcpp::cppFunction(depends = "RcppArmadillo", "
arma::sp_mat f_cpp(const arma::mat &x, const arma::vec &y, double alpha, Rcpp::Function f_R) {
    arma::sp_mat coef = Rcpp::as<arma::sp_mat>(f_R(x, y, alpha));
    return coef;
}")
coefs2 <- f_cpp(X, y, alpha = 1, my.glmnet)

all(coefs - coefs2 == 0)
#> [1] TRUE

Создано в 2019-06-12 пакетом Представления (v0.3.0)

Конечно, вы можете делать более интересные вещи с вычисленными коэффициентами, чем возвращать их в R. Явное Rcpp::as необходимо, поскольку в C ++ нет способа узнать, какой тип аргумента возвращает функция R. В этом случае это разреженная матрица, которую можно преобразовать в arma::sp_mat. Кстати, это приводит к потере Dimnames матрицы, поэтому нельзя использовать all.equal для сравнения.

...