R накопительная сумма с условием - PullRequest
0 голосов
/ 20 ноября 2018

(Для тех из вас, кто знаком с 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 Ответ

0 голосов
/ 21 ноября 2018

Более быстрая версия 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
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...