R中定积分的Monte-Carlo方法

Question

我开始了一个Msc,我们正在学习R包,但是我在运动时遇到了问题.练习是这样的:

"假设我们想要使用基本的蒙特卡罗方法估计∫x2(在1和0之间确定).基本上,我们在曲线上投掷飞镖并计算落在曲线下方的飞镖数量.该方法的算法包括:

1)初始化:命中= 0

2)for(i in 1:N):生成两个随机数,U1,U2在0和1之间.如果U2 <(U1)^ 2,则命中=命中+ 1结束.

3)面积估计=命中/ N.

使用循环提供R代码以实现此蒙特卡罗算法.你的功能需要多长时间？提供更有效的代码,避免以前的循环.说明通过矢量化代码可以提高效率."

我有这些代码,但我认为我做错了.

montecarlo <- function(N){
  hits=0
  for (i in 1:N){
    U1 <- runif(1, 0, 1)
    U2 <- runif(1, 0, 1)
    if (U2 < (U1)^2){
      hits = hits+1}}
  return(hits/N)
}

montecarlo2 <- function(N){
  hits=0
  U1 <- runif (1:N, 0, 1)
  U2 <- runif (1:N, 0, 1)
  hits= hits+1 [U2<(U1)^2]
  return(hits/N)
}

对于第一种方法,使用循环,我已经获得(例如):

> montecarlo(23)
[1] 0.3478261
> montecarlo(852)
[1] 0.3591549
> montecarlo(8563255)
[1] 0.3332472

你能帮助我吗？非常感谢你:S.

Answer 1

其中一种方式:

montecarlo_for <- function(N) {
  hits <- 0
  for (i in 1:N) {
    U1 <- runif(1)
    U2 <- runif(1)
    if (U2 < (U1) ^ 2) hits <- hits + 1
  }
  return(hits / N)
}

矢量化

montecarlo_vec <- function(N) {
  sum(runif(N) < runif(N)^2) / N
}

比较速度,例如使用microbenchmark包:

microbenchmark::microbenchmark(times = 50,
  montecarlo_for(1e5),
  montecarlo_vec(1e5)
)

我的机器上的速度比较表明,矢量化方法大约快了100倍(平均值和中位数时间如下所示):

Unit: milliseconds
expr                  mean       median
montecarlo_for(1e+05) 509.927001 497.238904
montecarlo_vec(1e+05) 5.214527   4.922007

只是为了好玩,如果你想看看算法收敛到结果(1/3)的速度越来越快,样本量越来越大:

plot(sapply(1:1000, montecarlo_vec), type = "line")