Rcpp 中的向量均值

Question

我正在尝试将 r 函数转换为 Rcpp 以尝试加快速度，因为它涉及 for 循环。一路上我需要计算向量条目的平均值，这在 R 中就像 mean(x) 一样简单，但它似乎在 Rcpp 中不起作用，每次都给我 0 0 作为结果。

我的代码如下所示：

cppFunction(
  "NumericVector fun(int n, double lambda, ...) {
   ...
   NumericVector y = rpois(n, lambda);
   NumericVector w = dpois(y, lambda);
   NumericVector x = w*y;
   double z = mean(x);
   return z;
}")

编辑：所以我认为我的错误是由于上面提到的，而 return 的一个 double z 只是我试图找出问题所在。然而，以下代码仍然不起作用：

cppFunction(
   "NumericVector zstat(int n, double lambda, double lambda0, int m) {
   NumericVector z(m);
   for (int i=1; i<m; ++i){
   NumericVector y = rpois(n, lambda0);
   NumericVector w = dpois(y, lambda)/dpois(y,lambda0);
   double x = mean(w*y);
   z[i] = (x-2)/(sqrt(2/n));
   }
   return z;
}")

Answer 1

您函数的 return 类型是 NumericVector，但 Rcpp::mean return 是可转换为 double 的标量值。修复此问题将解决问题：

library(Rcpp)

cppFunction(
  "double fun(int n, double lambda) {
   NumericVector y = rpois(n, lambda);
   NumericVector w = dpois(y, lambda);
   NumericVector x = w*y;
   double z = mean(x);
   return z;
}")

set.seed(123)
fun(50, 1.5)
# [1] 0.2992908

您的代码中发生的事情是因为 NumericVector 被指定为 return 类型，this constructor is called、

template <typename T>
Vector(T size, 
    typename Rcpp::traits::enable_if<traits::is_arithmetic<T>::value, void>::type* = 0) {
    Storage::set__( Rf_allocVector( RTYPE, size) ) ;
    init() ;
}

将 double 转换为整数类型并创建长度等于 double 的截断值的 NumericVector。为了演示，

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector from_double(double x) {
    return x;
}

/*** R

sapply(0.5:4.5, from_double)
# [[1]]
# numeric(0)
#
# [[2]]
# [1] 0
#
# [[3]]
# [1] 0 0
#
# [[4]]
# [1] 0 0 0
#
# [[5]]
# [1] 0 0 0 0

*/

编辑：关于你的第二个问题，你除以sqrt(2 / n)，其中2和n都是整数，结束在大多数情况下，向上导致除以零——因此结果向量中的所有 Inf 值。您可以使用 2.0 而不是 2 来解决此问题：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector zstat(int n, double lambda, double lambda0, int m) {
    NumericVector z(m);
    for (int i=1; i<m; ++i){
        NumericVector y = rpois(n, lambda0);
        NumericVector w = dpois(y, lambda)/dpois(y,lambda0);
        double x = mean(w * y);
        // z[i] = (x - 2) / sqrt(2 / n);
        //                       ^^^^^
        z[i] = (x - 2) / sqrt(2.0 / n);
        //                    ^^^^^^^
   }
   return z;
}

/*** R

set.seed(123)
zstat(25, 2, 3, 10)
# [1]  0.0000000 -0.4427721  0.3199805  0.1016661  0.4078687  0.4054078
# [7] -0.1591861  0.9717596  0.6325110  0.1269779

*/

C++ 不是 R -- 您需要更加注意变量的类型。

Rcpp 中的向量均值

Vector mean in Rcpp

r

rcpp