【问题标题】:Why do I get the error for using "pnorm" in Rcpp为什么在 Rcpp 中使用“pnorm”时出现错误
【发布时间】:2018-11-14 20:50:23
【问题描述】:

我需要在我的 Rcpp 代码中包含来自 arma:: 的变量。但是我在尝试使用糖函数pnorm 时遇到了问题。这是一个演示:

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

// [[Rcpp::export]]
double pget(NumericVector x, NumericVector beta) {
  arma::colvec xx = Rcpp::as<arma::colvec>(x) ;
  arma::colvec bb = Rcpp::as<arma::colvec>(beta) ;
  double tt = as_scalar( arma::trans(xx) * bb);
  double temp = Rcpp::pnorm(tt);
  return temp;
}

然后报错:no matching function for call to 'pnorm5'

这是否意味着我不能使用Rcpp::pnorm???

【问题讨论】:

  • 这几乎肯定是重复了几次,但我现在没有时间去挖掘。

标签: r rcpp armadillo


【解决方案1】:

Rcpp 糖函数适用于矢量类型参数,例如Rcpp::NumericVector。对于标量参数,您可以使用 R 命名空间中的函数:

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

// [[Rcpp::export]]
double pget(NumericVector x, NumericVector beta) {
  arma::colvec xx = Rcpp::as<arma::colvec>(x) ;
  arma::colvec bb = Rcpp::as<arma::colvec>(beta) ;
  double tt = as_scalar( arma::trans(xx) * bb);
  double temp = R::pnorm(tt, 0.0, 1.0, 1, 0);
  return temp;
}

/*** R
x <- rnorm(5)
beta <- rnorm(5)
pget(x, beta)
*/

顺便说一句,这里有两种变体。第一个变体使用arma 而不是Rcpp 向量作为参数。由于这些是const 引用,因此不会复制任何数据。另外,使用arma::dot

// [[Rcpp::export]]
double pget2(const arma::colvec& xx, const arma::colvec& bb) {
  double tt = arma::dot(xx, bb);
  return R::pnorm(tt, 0.0, 1.0, 1, 0);
}

第二个变体计算标量积而不使用犰狳:

// [[Rcpp::export]]
double pget3(NumericVector x, NumericVector beta) {
  double tt = Rcpp::sum(x * beta);
  return R::pnorm(tt, 0.0, 1.0, 1, 0);
}

【讨论】:

  • 谢谢!如果我进行矢量化,Rcpp::pnorm 会比R::pnorm 快吗?
  • 就解决 OP 的问题而言,答案很好 (+1),但请注意您发布的代码:您的参数顺序有些混乱。你想要R::pnorm(tt, 0.0, 1.0, 1, 0)而不是R::pnorm(tt, 1.0, 0.0, 1, 0)(平均值在sd之前)。
  • @lostintheforest 对这两个函数进行基准测试怎么样?
  • as_scalar(arma::trans(xx) * bb) 不等于更短的 dot(xx,bb) 吗? (documentation)
  • @mtall 确实,感谢您的提醒。我已经更新了我的答案。
【解决方案2】:

我比 Rcpp 的 @RalfStubner 更不是专家,所以我不得不四处寻找(在 StackOverflowRcpp cheat sheat 的帮助下)来获取以下代码。我没有在标量上使用 R 命名空间版本,而是转换回 NumericVector ......这几乎可以肯定地由真正知道自己在做什么的人更有效地完成/跳过几步......例如有可能不经过as_scalar 就可以直接完成 arma 到 NumericVector 的转换...?

#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>
#include <Rcpp.h>

// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
NumericVector pget(NumericVector x, NumericVector beta) {
  colvec xx = as<colvec>(x) ;
  colvec bb = as<colvec>(beta) ;
  double tt = as_scalar(trans(xx) * bb);
  NumericVector tt2 = NumericVector::create( tt );
  NumericVector temp = Rcpp::pnorm(tt2);
  return temp;
}

【讨论】:

  • armaNumericVector 的转换可以在不经过 as_scalar() 的情况下完成,如下所示:NumericVector tt = as&lt;NumericVector&gt;(wrap(arma::trans(xx) * bb));,但我还没有进行基准测试看看它是否有效
猜你喜欢
  • 1970-01-01
  • 2016-07-30
  • 1970-01-01
  • 2021-06-26
  • 1970-01-01
  • 2021-04-18
  • 2021-04-20
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多