# \donttest{
# This example shows how to use the optimizers
# for C++ objective functions. We will use
# a linear regression as an example. Note that
# this is not a useful application of the optimizers
# as there are specialized packages for linear regression
# (e.g., glmnet)
library(Rcpp)
library(lessSEM)
linreg <- '
// [[Rcpp::depends(RcppArmadillo)]]
#include
// [[Rcpp::export]]
double fitfunction(const Rcpp::NumericVector& parameters, Rcpp::List& data){
// extract all required elements:
arma::colvec b = Rcpp::as(parameters);
arma::colvec y = Rcpp::as(data["y"]); // the dependent variable
arma::mat X = Rcpp::as(data["X"]); // the design matrix
// compute the sum of squared errors:
arma::mat sse = arma::trans(y-X*b)*(y-X*b);
// other packages, such as glmnet, scale the sse with
// 1/(2*N), where N is the sample size. We will do that here as well
sse *= 1.0/(2.0 * y.n_elem);
// note: We must return a double, but the sse is a matrix
// To get a double, just return the single value that is in
// this matrix:
return(sse(0,0));
}
// [[Rcpp::export]]
arma::rowvec gradientfunction(const Rcpp::NumericVector& parameters, Rcpp::List& data){
// extract all required elements:
arma::colvec b = Rcpp::as(parameters);
arma::colvec y = Rcpp::as(data["y"]); // the dependent variable
arma::mat X = Rcpp::as(data["X"]); // the design matrix
// note: we want to return our gradients as row-vector; therefore,
// we have to transpose the resulting column-vector:
arma::rowvec gradients = arma::trans(-2.0*X.t() * y + 2.0*X.t()*X*b);
// other packages, such as glmnet, scale the sse with
// 1/(2*N), where N is the sample size. We will do that here as well
gradients *= (.5/y.n_rows);
return(gradients);
}
// Dirk Eddelbuettel at
// https://gallery.rcpp.org/articles/passing-cpp-function-pointers/
typedef double (*fitFunPtr)(const Rcpp::NumericVector&, //parameters
Rcpp::List& //additional elements
);
typedef Rcpp::XPtr fitFunPtr_t;
typedef arma::rowvec (*gradientFunPtr)(const Rcpp::NumericVector&, //parameters
Rcpp::List& //additional elements
);
typedef Rcpp::XPtr gradientFunPtr_t;
// [[Rcpp::export]]
fitFunPtr_t fitfunPtr() {
return(fitFunPtr_t(new fitFunPtr(&fitfunction)));
}
// [[Rcpp::export]]
gradientFunPtr_t gradfunPtr() {
return(gradientFunPtr_t(new gradientFunPtr(&gradientfunction)));
}
'
Rcpp::sourceCpp(code = linreg)
ffp <- fitfunPtr()
gfp <- gradfunPtr()
N <- 100 # number of persons
p <- 10 # number of predictors
X <- matrix(rnorm(N*p), nrow = N, ncol = p) # design matrix
b <- c(rep(1,4),
rep(0,6)) # true regression weights
y <- X%*%matrix(b,ncol = 1) + rnorm(N,0,.2)
data <- list("y" = y,
"X" = cbind(1,X))
parameters <- rep(0, ncol(data$X))
names(parameters) <- paste0("b", 0:(length(parameters)-1))
al1 <- gpAdaptiveLassoCpp(par = parameters,
regularized = paste0("b", 1:(length(b)-1)),
fn = ffp,
gr = gfp,
lambdas = seq(0,1,.1),
additionalArguments = data)
al1@parameters
# }
Run the code above in your browser using DataLab