使用Rcpp在C++中的R中应用optim函数

Question

我试图在Rcpp中调用R函数.我在C++中使用Rcpp看到了一个调用R的optim函数的例子,但是我无法正确地修改它用于我的用例.基本上,目标函数取决于并且我想要优化它.optim()xyb

这是执行我想要的R代码:

example_r = function(b, x, y) {
  phi = rnorm(length(x))

  tar_val = (x ^ 2 + y ^ 2) * b * phi

  objftn_r = function(beta, x, y) {
    obj_val = (x ^ 2 + y ^ 2) * beta

    return(obj_val)
  }

  b1 = optim(b, function(beta) {
    sum((objftn_r(beta, x, y) - tar_val) ^ 2)
  }, method = "BFGS")$par

  result = (x ^ 2 + y ^ 2) * b1

  return(b1)
}

这是我尝试将其翻译为_RcppArmadillo:

#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::depends(RcppArmadillo)]]

arma::vec example_rcpp(arma::vec b, arma::vec x, arma::vec y){

  arma::vec tar_val = pow(x,2)%b-pow(y,2);

  return tar_val;
}

// [[Rcpp::export]]
arma::vec optim_rcpp(const arma::vec& init_val, arma::vec& x, arma::vec& y){

  Rcpp::Environment stats("package:stats"); 
  Rcpp::Function optim = stats["optim"];

  Rcpp::List opt_results = optim(Rcpp::_["par"]    = init_val,
                                 Rcpp::_["fn"]     = Rcpp::InternalFunction(&example_rcpp),
                                 Rcpp::_["method"] = "BFGS");

  arma::vec out = Rcpp::as<arma::vec>(opt_results[0]);

  return out;
} 

但是,此代码返回:

> optim_rcpp(1:3,2:4,3:5)
Error in optim_rcpp(1:3, 2:4, 3:5) : not compatible with requested type

我不确定这里的错误是什么.

Answer 1

在开始之前,我有几点意见:

1. 请显示您的所有尝试.
   • 特别是,请确保您的示例是一个可重现性最小的示例
2. 不要不删除或缩短代码,除非问.
3. 保持问题的范围狭窄.
   • 使用optim从- [R在C++比使用以非常不同的C++底层C++代码opt()从nlopt.
4. 避免垃圾邮件问题.
   • 如果您发现自己快速连续询问了3个以上的问题,请阅读文档或亲自与熟悉内容的人员交谈.

因此,我已经清理了你的问题......但是,将来,这种情况很可能不会发生.

数据生成过程

数据生成过程似乎分两步完成:首先,在example_r函数外部,然后在函数内部.

这应该简化,以便在优化功能之外完成.例如:

generate_data = function(n, x_mu = 0, y_mu = 1, beta = 1.5) {

  x = rnorm(n, x_mu)
  y = rnorm(n, y_mu)

  phi = rnorm(length(x))

  tar_val = (x ^ 2 + y ^ 2) * beta * phi

  simulated_data = list(x = x, y = y, beta = beta, tar_val = tar_val)
  return(simulated_data)
}

客观函数和R 'optim

目标函数必须在R中返回单个值,例如标量.在已发布的R代码下,实际上有两个功能被设计为按顺序充当目标函数,例如

objftn_r = function(beta, x, y) {
  obj_val = (x ^ 2 + y ^ 2) * beta

  return(obj_val)
}

b1 = optim(b, function(beta) {
  sum((objftn_r(beta, x, y) - tar_val) ^ 2)
}, method = "BFGS")$par

因此,该目标函数应重写为:

objftn_r = function(beta_hat, x, y, tar_val) {

  # The predictions generate will be a vector
  est_val = (x ^ 2 + y ^ 2) * beta_hat

  # Here we apply sum of squares which changes it
  # from a vector into a single "objective" value
  # that optim can work with.
  obj_val = sum( ( est_val  - tar_val) ^ 2)

  return(obj_val)
}

从那里,调用应该对齐为:

sim_data = generate_data(10, 1, 2, .3)

b1 = optim(sim_data$beta, fn = objftn_r, method = "BFGS",
           x = sim_data$x, y = sim_data$y, tar_val = sim_data$tar_val)$par

RcppArmadillo目标函数

修正了R代码的范围和行为后,让我们专注于将其转换为RcppArmadillo.

特别是,注意到,翻译后所定义的目标函数返回一个向量而不是标量到optim,这是不一个单一的值.另外值得关注的是tar_val目标函数中缺少参数.考虑到这一点,目标函数将转化为:

// changed function return type and 
// the return type of first parameter
double obj_fun_rcpp(double& beta_hat, 
                    arma::vec& x, arma::vec& y, arma::vec& tar_val){

  // Changed from % to * as it is only appropriate if  
  // `beta_hat` is the same length as x and y.
  // This is because it performs element-wise multiplication
  // instead of a scalar multiplication on a vector
  arma::vec est_val = (pow(x, 2) - pow(y, 2)) * beta_hat;

  // Compute objective value
  double obj_val = sum( pow( est_val - tar_val, 2) );

  // Return a single value
  return obj_val;
}

现在,目标函数集,让我们解决RCPP调用放入[R为optim()从C++.在此函数中,必须显式提供函数的参数.所以x,y和tar_val必须存在的optim呼叫.因此,我们将最终得到:

// [[Rcpp::export]]
arma::vec optim_rcpp(double& init_val,
                     arma::vec& x, arma::vec& y, arma::vec& tar_val){

  // Extract R's optim function
  Rcpp::Environment stats("package:stats"); 
  Rcpp::Function optim = stats["optim"];

  // Call the optim function from R in C++ 
  Rcpp::List opt_results = optim(Rcpp::_["par"]    = init_val,
                                 // Make sure this function is not exported!
                                 Rcpp::_["fn"]     = Rcpp::InternalFunction(&obj_fun_rcpp),
                                 Rcpp::_["method"] = "BFGS",
                                 // Pass in the other parameters as everything
                                 // is scoped environmentally
                                 Rcpp::_["x"] = x,
                                 Rcpp::_["y"] = y,
                                 Rcpp::_["tar_val"] = tar_val);

  // Extract out the estimated parameter values
  arma::vec out = Rcpp::as<arma::vec>(opt_results[0]);

  // Return estimated values
  return out;
}

全部一起

完整功能的代码可以通过以下方式编写test_optim.cpp和编译sourceCpp():

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

// changed function return type and 
// the return type of first parameter
// DO NOT EXPORT THIS FUNCTION VIA RCPP ATTRIBUTES
double obj_fun_rcpp(double& beta_hat, 
                    arma::vec& x, arma::vec& y, arma::vec& tar_val){

  // Changed from % to * as it is only appropriate if  
  // `beta_hat` is the same length as x and y.
  // This is because it performs element-wise multiplication
  // instead of a scalar multiplication on a vector
  arma::vec est_val = (pow(x, 2) - pow(y, 2)) * beta_hat;

  // Compute objective value
  double obj_val = sum( pow( est_val - tar_val, 2) );

  // Return a single value
  return obj_val;
}


// [[Rcpp::export]]
arma::vec optim_rcpp(double& init_val,
                     arma::vec& x, arma::vec& y, arma::vec& tar_val){

  // Extract R's optim function
  Rcpp::Environment stats("package:stats"); 
  Rcpp::Function optim = stats["optim"];

  // Call the optim function from R in C++ 
  Rcpp::List opt_results = optim(Rcpp::_["par"]    = init_val,
                                 // Make sure this function is not exported!
                                 Rcpp::_["fn"]     = Rcpp::InternalFunction(&obj_fun_rcpp),
                                 Rcpp::_["method"] = "BFGS",
                                 // Pass in the other parameters as everything
                                 // is scoped environmentally
                                 Rcpp::_["x"] = x,
                                 Rcpp::_["y"] = y,
                                 Rcpp::_["tar_val"] = tar_val);

  // Extract out the estimated parameter values
  arma::vec out = Rcpp::as<arma::vec>(opt_results[0]);

  // Return estimated values
  return out;
}

测试用例

# Setup some values
beta = 2
x = 2:4
y = 3:5

# Set a seed for reproducibility
set.seed(111)

phi = rnorm(length(x))

tar_val = (x ^ 2 + y ^ 2) * beta * phi

optim_rcpp(beta, x, y, tar_val)
#          [,1]
# [1,] 2.033273

注意:如果您想避免返回大小为1 x1的矩阵,请使用double作为返回参数optim_rcpp并切换Rcpp::as<arma::vec>到Rcpp::as<double>