问题描述
我正在尝试使用简单的递归实现来计算变量的平均值:
I am trying to compute the mean of a variable using a simple, recursive implementation:
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 中实现类似的东西:
Here, m
will quickly converge to 2, which corresponds to the mean of the normal distribution where we generate new data points from. Implementing something similar in 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 的翻译错误在哪里吗?
This code does not show the expected behavior. Instead, it returns the initial value. Can somebody enlighten me where my error in this translation from R to Rcpp is?
推荐答案
在这一行:
m = ((irep)/(irep+1)) * m + (1/(irep+1)) * new_data;
您将一个 int
除以另一个 int
两次.在 C++ 中,整数除法返回另一个整数,丢弃余数.
You are dividing an int
by another int
twice. In C++, integer division returns another integer, discarding the remainder.
为了得到你想要的,强制以浮点形式进行除法:
To get what you want, force the divisions to be done in floating point:
m = ((irep)/(irep+1.0)) * m + (1.0/(irep+1.0)) * new_data;
这篇关于R 与 Rcpp 中的递归均值的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!