(对于那些熟悉 MCMC 的人,我正在尝试编写 Metropolis-Hastings 算法的(一个步骤))。
我正在尝试计算起始值为 0.5 的小随机值向量的累加和。但是,如果任何一点的累计总和小于 0 或大于 1,我需要复制之前的值并继续累计总和,而不对这些值求和,这将打破这种情况。
注意:我需要一个向量化的解决方案(无循环或索引)以进行优化或快速处理。仅使用基本 R 函数的奖励积分。
例子:
set.seed(1)
temp=c(0.5,runif(20,-0.3,0.3))
cumsum(temp)
[1] 0.5000000 0.3593052 0.2825795 0.3262916 0.5712162 0.3922254 0.6312592
[8] 0.8980644 0.9945430 1.0720115 0.8090832 0.6326680 0.4386020 0.5508157
[15] 0.4812780 0.6431828 0.6418024 0.7723735 1.0675171 0.9955382 1.1620054
但是我需要的是
[1] 0.5000000 0.3593052 0.2825795 0.3262916 0.5712162 0.3922254 0.6312592
[8] 0.8980644 0.9945430 0.9945430 0.7316148 0.5551995 0.3611336 0.4733473
[15] 0.4038095 0.5657144 0.5643339 0.6949050 0.9900487 0.9180698 0.9180698
我们可以使用 for 循环来做到这一点
for (i in 2:21) {
temp[i]=temp[i-1]+temp[i]
if(temp[i]<0 | temp[i]>1) {
temp[i]=temp[i-1]
}
}
最佳答案
更快的 C++ 版本:
library(Rcpp)
Cpp_boundedCumsum <- cppFunction('NumericVector boundedCumsum(NumericVector x){
int n = x.size();
NumericVector out(n);
double tmp;
out[0] = x[0];
for(int i = 1; i < n; ++i){
tmp = out[i-1] + x[i];
if(tmp < 0.0 || tmp > 1.0)
out[i] = out[i-1];
else
out[i] = tmp;
}
return out;
}')
与R版本对比:
R_boundedCumsum <- function(x){
for (i in 2:length(x)){
x[i] <- x[i-1]+x[i]
if(x[i]<0 || x[i]>1)
x[i] <- x[i-1]
}
x
}
x <- runif(1000)
all.equal(R_boundedCumsum(x), Cpp_boundedCumsum(x))
[1] TRUE
library(microbenchmark)
microbenchmark(R_boundedCumsum(x), Cpp_boundedCumsum(x))
Unit: microseconds
expr min lq mean median uq max neval
R_boundedCumsum(x) 2062.629 2262.2225 2460.65661 2319.358 2562.256 4112.540 100
Cpp_boundedCumsum(x) 3.636 4.3475 7.06454 5.792 9.127 25.703 100
关于带条件的 R 累积和,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/53398554/