Why do I get the error for using "pnorm" in Rcpp

Question

I need to involve variable from arma::in my Rcpp code. But I ran into a problem when trying to use the sugar function pnorm. Here is a demo:

#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;
}

Then I got an error: no matching function for call to 'pnorm5'

Does that mean I cannot use Rcpp::pnorm???

This almost surely is a duplicate a few times over but I don't have time to dig now. — Dirk Eddelbuettel, Nov 14 '18 at 21:21

Ralf Stubner · Accepted Answer · 2018-11-15T10:00:20.967

7

The Rcpp sugar functions are meant for vector type arguments like Rcpp::NumericVector. For scalar arguments you can use the functions in the R namespace:

#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)
*/

BTW, here two variants. First variant uses arma instead of Rcpp vectors as arguments. Since these are const references, no data is copied. In addition, arma::dot is used:

// [[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);
}

The second variant calculates the scalar product without resorting to Armadillo:

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

edited Nov 15 '18 at 10:00

answered Nov 14 '18 at 21:04

Ralf Stubner

26,263
3
40
75

Thanks! Will `Rcpp::pnorm` be faster than `R::pnorm` if I do vectorization? – lostintheforest Nov 14 '18 at 21:08
2

Great answer (+1) in terms of getting down to OP's problem, but a note on your posted code: You have the order of arguments mixed up a little. You want `R::pnorm(tt, 0.0, 1.0, 1, 0)` rather than `R::pnorm(tt, 1.0, 0.0, 1, 0)` (mean comes before sd). – duckmayr Nov 14 '18 at 21:28
@lostintheforest How about benchmarking the two functions? – Ralf Stubner Nov 14 '18 at 23:27
2

isn't `as_scalar(arma::trans(xx) * bb)` equivalent to the much shorter `dot(xx,bb)` ? ([documentation](http://arma.sourceforge.net/docs.html#dot)) – mtall Nov 15 '18 at 01:20
@mtall Indeed, thanks for the reminder. I have updated my answer. – Ralf Stubner Nov 15 '18 at 10:01

score 2 · Answer 2 · answered Nov 14 '18 at 21:25

I'm much less of an expert than @RalfStubner at Rcpp, so I had to hack around (with help from StackOverflow and the Rcpp cheat sheat) to get the following code. Instead of using the R-namespace versions on scalars, I converted back to a NumericVector ... this can almost certainly be done more efficiently/skipping a few steps by someone who actually knows what they're doing ... e.g. it's possible that the arma-to-NumericVector conversion could be done directly without going through as_scalar ... ?

#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;
}

The `arma` to `NumericVector` conversion can be done without going through `as_scalar()`, like this: `NumericVector tt = as(wrap(arma::trans(xx) * bb));`, but I haven't benchmarked to see if it's efficient — duckmayr, Nov 14 '18 at 21:38

Why do I get the error for using "pnorm" in Rcpp

2 Answers2