代码之家  ›  专栏  ›  技术社区  ›  user2974951

r有条件累积和

  •  2
  • user2974951  · 技术社区  · 6 年前

    (对于你们中熟悉MCMC的人,我正在尝试编写大都会黑斯廷斯算法的一个步骤)。

    我试图做一个小随机值向量的累积和,初始值为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]
        }
    }
    
    1 回复  |  直到 6 年前
        1
  •  8
  •   Nairolf    6 年前

    一个更快的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