Rcppでループをまわす話

inspired by http://bicycle1885.hatenablog.com/entry/2015/01/04/175916

速度がどうというより、愚直なループをみたらRcppをさくっと適用できるようになりたい。

見本(Julia)はこちら。

circle_in = 0.0
for i in 1:100000000
  l = (rand()^2 + rand()^2) ^ 0.5
  if l <= 1
    circle_in = circle_in + 1
  end
end
println(4 * circle_in / 100000000)

Rで愚直に書くとこんな感じ。

calcPiR <- function(){
  circle_in = 0
  for (i in 1:100000000){
    l = sqrt(runif(1)^2 + (runif(1))^2)
    if (l <= 1){
      circle_in = circle_in + 1      
    }
  }
  4 * circle_in / 100000000 
}

しかしこれだといつまで経っても終わらないので少しはましなコードにしてそちらを実験に使う。

ましなコードはこちらからもってきた。

http://gallery.rcpp.org/articles/simulating-pi/

ということでRと愚直なC++のループとRのコードっぽいC++で比較した。

library(Rcpp)
library(microbenchmark)

calcPiR2 <- function() {
    x <- runif(100000000)
    y <- runif(100000000)
    d <- sqrt(x^2 + y^2)
    return(4 * sum(d < 1.0) / 100000000)
}

cppFunction('double calcPiCpp() {
            double circle_in = 0.0;
              for (int i = 0; i < 100000000; i++) {
                double l = sqrt(pow(rand()/2147483648., 2)+pow(rand()/2147483648.,2));
                if (l <= 1){
                  circle_in = circle_in + 1.0;
                  }
              }
              return 4*circle_in/100000000;
            }')

cppFunction('double calcPiCpp2() {
            NumericVector x = runif(100000000);
            NumericVector y = runif(100000000);
            NumericVector d = sqrt(x*x + y*y);
            return 4.0 * sum(d < 1.0) / 100000000;
            }')

# @kos59125より
calc_pi <- function(n, k) {
   stopifnot(n %% k == 0)
   r <- n %/% k
   s <- 0
   loop <- 0L
   while (loop < r) {
      x <- runif(k)
      y <- runif(k)
      s <- s + sum(x * x + y * y <= 1)
      loop <- loop + 1L
   }
   4 * s / n
}
 
microbenchmark(calcPiR2(), calcPiCpp(), calcPiCpp2(), calc_pi(100000000L, 10000L), times=3L)

## Unit: seconds
##                         expr       min        lq      mean    median
##                   calcPiR2() 10.709923 10.710115 11.513909 10.710307
##                  calcPiCpp()  2.355628  2.358802  2.374402  2.361975
##                 calcPiCpp2()  5.680846  5.695607  5.791551  5.710368
##  calc_pi(100000000L, 10000L)  9.068039  9.424003  9.547543  9.779967
##         uq       max neval
##  11.915902 13.121498     3
##   2.383789  2.405603     3
##   5.846904  5.983440     3
##   9.787296  9.794624     3

環境はこちら。

sessionInfo()

## R version 3.1.1 (2014-07-10)
## Platform: x86_64-apple-darwin13.1.0 (64-bit)
## 
## locale:
## [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
## 
## attached base packages:
## [1] stats     grDevices datasets  graphics  utils     methods   base     
## 
## other attached packages:
## [1] microbenchmark_1.4-2 Rcpp_0.11.3         
## 
## loaded via a namespace (and not attached):
##  [1] audio_0.1-5      beepr_1.1        colorspace_1.2-4 digest_0.6.4    
##  [5] evaluate_0.5.5   formatR_1.0      ggplot2_1.0.0    grid_3.1.1      
##  [9] gtable_0.1.2     htmltools_0.2.6  knitr_1.7        MASS_7.3-33     
## [13] munsell_0.4.2    plyr_1.8.1       proto_0.3-10     reshape2_1.4    
## [17] rmarkdown_0.3.10 scales_0.2.4     stringr_0.6.2    tools_3.1.1     
## [21] yeah_0.3.3