R 与 Rcpp 中的递归均值

Question

我正在尝试使用简单的递归实现来计算变量的平均值：

m <- 0 # initialize mean

for(irep in 0:999){

# new data point
new_data <- rnorm(1,2,1)

# recursive formula for sample mean
m = (irep/(irep+1)) * m + (1/(irep+1)) * new_data

}

在这里，m将快速收敛到 2，这对应于我们从中生成新数据点的正态分布的均值。在 Rcpp 中实现类似的东西：

#include <RcppArmadillo.h>

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

using namespace Rcpp;

// [[Rcpp::export]]
double  rec_mean(int sample_size){

  double m = 0; //initialize mean

  for(int irep = 0; irep < sample_size; irep++){

      // new data
      double new_data = R::rnorm(2,1);

      // mean recursive update
      m = ((irep)/(irep+1)) * m + (1/(irep+1)) * new_data;

  }

  return m;

}

此代码未显示预期行为。相反，它返回初始值。有人能告诉我我从 R 到 Rcpp 的翻译错误在哪里吗？

Answer 1

在这一行：

m = ((irep)/(irep+1)) * m + (1/(irep+1)) * new_data;

您将一个int除以另一个int两次。在 C++ 中，整数除法返回另一个整数，丢弃余数。

为了得到你想要的，强制除法以浮点形式完成：

m = ((irep)/(irep+1.0)) * m + (1.0/(irep+1.0)) * new_data;