R中定积分的蒙特卡洛方法

Question

我开始了一个 Msc，在那里我们正在学习 R 包，但我在练习时遇到了问题。练习是这样的：

"假设我们要使用基本的蒙特卡洛方法估计 ∫x2（确定在 1 和 0 之间）。本质上，我们在曲线上扔飞镖并计算落在曲线下方的飞镖的数量。算法该方法包括：

1) 初始化：hits=0

2) for (i in 1:N): 生成两个0和1之间的随机数，U1,U2。如果U2<(U1)^2，则hits=hits+1 end for.

3) 面积估计 = hits/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")

Answer 2

我发现这是一种有助于估算积分的方法。考虑以下因素：

$\int_A f(x)dx = \int_A \frac{f(x)}{g(x)}g(x)dx = E_g\left(\frac{f(x)}{g(x)}\right)$

所以如果我们在支持 A 上有一个分布 G，那么我们可以用分布 G 中的样本估计 f(x) 在 A 上的积分。函数 f 不必严格为正。

在你的例子中，让X均匀分布，那么

$\int_0^1 x^2dx = \int_0^1 \frac{x^2}{g(x)}g(x)dx = E_g\left(\frac{x^2}{g(x)}\right) = E_g(x^2)$

因为对于支持中的所有 x，g(x)=1。所以你可以估计积分

N = 100000
mean(runif(N)^2)

你也可以让 X 成为一个 beta 随机变量并用

估计积分

x = rbeta(N, 2, 1)
fx = x^2
gx = dbeta(x, 2, 1)
mean(fx / gx)

这是另一个在正实数线上估计函数的例子。

f = function(x)
    abs(sin(x)) / (x+1)^2
x = rgamma(N, 2, 1/4)
fx = f(x)
gx = dgamma(x, 2, 1/4)
mean(fx / gx)

我不记得所有的理论（对于何时可以使用这种方法可能会有一些限制），但通常如果你可以在与积分相同的支持下从任何分布中抽样，你可以估计积分如上。如果你选择的 G 对函数 f 的密度是 "closer"，那么你的 Monte Carlo 误差就会越小。

R中定积分的蒙特卡洛方法

Monte-Carlo method for definite integral in R

r

montecarlo