I am trying to make some statistical calculations in R faster by using Rcpp. First, I wrote the code in R. Then I wrote the code in C using Qt Creator, which required me to use the Boost package. Now, when I try to use sourceCpp() to compile a simple function that uses boost::math::statistics::two_sample_t_test(), I get two errors--click here to see how it looks in RStudio.

/Library/Frameworks/R.framework/Versions/3.6/Resources/library/BH/include/boost/compute/algorithm/random_shuffle.hpp
no member named 'random_shuffle' in namespace 'std'; did you mean simply 'random_shuffle'?

~/Documents/Research/P-value correction project/Perm FDR with C using Rcpp/PermFDR_R/Rcpp.cpp
no member named 'two_sample_t_test' in namespace 'boost::math::statistics'

Here is a screenshot of the code in /Library/Frameworks/R.framework/Versions/3.6/Resources/library/BH/include/boost/compute/algorithm/random_shuffle.hpp that pops up when this happens.

Here is the R code.

library(Rcpp)
library(rstudioapi)
library(BH)

Sys.setenv("PKG_CXXFLAGS"="-std=c  17")
sourceCpp("Rcpp.cpp")

Here is the C code.

//[[Rcpp::depends(BH)]]
#include <Rcpp.h>
#include <vector>
#include <cmath>
#include <iostream>
#include <random>
#include <boost/math/statistics/t_test.hpp>
#include <boost/math/distributions/students_t.hpp>
#include <boost/math/tools/univariate_statistics.hpp>
#include <boost/compute/algorithm/random_shuffle.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <iomanip>
#include <numeric>
#include <random>

//[[Rcpp::plugins(cpp17)]]

using namespace std;
#include <Rcpp.h>
using namespace Rcpp;

/*
 * Performs a T test on the measurements according to the design (1s and 2s)
 * and returns a P value.
 */
double designTTest(vector<double> ints, vector<int> design) {
  if (ints.size() != design.size()) {
    cout << "ERROR: DESIGN VECTOR AND MEASUREMENT VECTOR NOT EQUAL IN LENGTH!";
    throw;
  }
  
  vector<double> cIntensities;
  vector<double> tIntensities;
  
  for (int i = 0; i < design.size(); i  ) {
    if (design[i] == 1) {
      cIntensities.push_back(ints[i]);
    } else if (design[i] == 2) {
      tIntensities.push_back(ints[i]);
    } else {
      cout << "ERROR: DESIGN SYMBOL IS NOT 1 OR 2!";
      throw;
    }
  }
  
  auto [t, p] = boost::math::statistics::two_sample_t_test(cIntensities, tIntensities);
  return p;
}

// [[Rcpp::export]]
double ttestC(NumericVector ints, NumericVector design) {
  vector<double> intVec = as<std::vector<double>>(ints);
  vector<int> designVec = as<std::vector<int>>(design);
  return designTTest(intVec, designVec);
}

Thank you!!

CodePudding user response：

There is a lot going on in your question, and I am not sure I understand all of (the 'design sorting' is very unclear).

I can, however, help you with the mechanics of Rcpp, and use of Boost via BH. I can suggest a number of changes:

• you do not need to specify either C 11 or C 17; R already defaults to C 14 (if the compiler supports it) under recent version
• you do not need all those headers: we need one for Rcpp, and one for Boost
• I recommend against using namespace ... and suggest explicit naming
• you do not need a wrapper from R vectors to std::vector<...> as Rcpp does that for you
• I simplified the error reporting via Rcpp::stop()

With all that, plus a mini-demo to run the function, your code becomes simpler and short.

Code

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

#include <Rcpp.h>
#include <boost/math/statistics/t_test.hpp>

// Performs a T test on the measurements according to the design (1s and 2s)
// and returns a P value.
// [[Rcpp::export]]
double designTTest(std::vector<double> ints, std::vector<double> design) {
    if (ints.size() != design.size()) Rcpp::stop("ERROR: DESIGN VECTOR AND MEASUREMENT VECTOR NOT EQUAL IN LENGTH!");
    std::vector<double> cIntensities, tIntensities;
    for (size_t i = 0; i < design.size(); i  ) {
        if (design[i] == 1) {
            cIntensities.push_back(ints[i]);
        } else if (design[i] == 2) {
            tIntensities.push_back(ints[i]);
        } else {
            Rcpp::stop("ERROR: DESIGN SYMBOL IS NOT 1 OR 2!");
        }
    }

    auto [t, p] = boost::math::statistics::two_sample_t_test(cIntensities, tIntensities);
    return p;
}

/*** R
designTTest(c(1,2,1,2,1), c(2,1,2,1,2))
*/

We can source this for compilation and the example (with possibly non-sensical data).

Output

> Rcpp::sourceCpp("~/git/stackoverflow/72652032/answer.cpp")

> designTTest(c(1,2,1,2,1), c(2,1,2,1,2))
[1] 0
> 

I removed a bunch of compilation noise that goes away when you add -Wno-parentheses to your CXXFLAGS in ~/.R/Makevars. CRAN does not let me (as BH maintainer) keep th upstream #pragmas so noisy it is ...

CodePudding user response：

In the full compiler output boost is telling you that you are not supposed to include <boost/math/tools/univariate_statistics.hpp>. You are supposed to include <boost/math/statistics/univariate_statistics.hpp> instead:

/usr/include/boost/math/tools/univariate_statistics.hpp:15:1: note: ‘#pragma message: This header is deprecated. Use <boost/math/statistics/univariate_statistics.hpp> instead.’

Following that advice, I don't get the std::random_shuffle error message anymore.

The other error message regarding two_sample_t_test is probably because that was added only relatively recently to boost, with version 1.76, see https://github.com/boostorg/math/pull/487. You might be using an older version.