我正在尝试使用简单的递归实现来计算变量的平均值:
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 的翻译错误在哪里吗?
最佳答案
在这一行中:
m = ((irep)/(irep+1)) * m + (1/(irep+1)) * new_data;
您将一个 int
除以另一个 int
两次。在 C++ 中,整数除法返回另一个整数,并丢弃余数。
为了得到你想要的,强制除法以 float 完成:
m = ((irep)/(irep+1.0)) * m + (1.0/(irep+1.0)) * new_data;
关于R 与 Rcpp 中的递归均值,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/61986077/