更高效的蒙特卡罗模拟制作

Question

所以，我已经编写了这段代码，可以有效地估计定义为h(x)的函数曲线下面积。我的问题是，我需要能够将面积估计到小数点后6位，但是我在estimateN中定义的算法似乎对我的计算机负担过重。基本上问题是如何使以下代码更加高效？有没有办法去掉循环？

h = function(x) {
    return(1+(x^9)+(x^3))
}
estimateN = function(n) {
    count = 0
    k = 1
    xpoints = runif(n, 0, 1)
    ypoints = runif(n, 0, 3)
    while(k <= n){
    if(ypoints[k]<=h(xpoints[k]))
        count = count+1
    k = k+1
    }
    #because of the range that im using for y
    return(3*(count/n))
}
#uses the fact that err<=1/sqrt(n) to determine size of dataset
estimate_to = function(i) {
    n = (10^i)^2
    print(paste(n, " repetitions: ", estimateN(n)))
}

estimate_to(6)

- Khodeir

---

estimate_to(6) 可能有点贪心：您当前的算法很可能会使 R 尝试分配长度为 1e12 的数值向量而耗尽内存。 - flodel

如果你真的需要进行1e12次模拟，那么你必须重写你的算法，以便在计算时间和内存使用之间找到一个折衷方案。 - flodel

使蒙特卡罗积分（更）高效的最佳方法，即在较少迭代次数内获得相同精度，是使用重要性采样，例如Metropolis Monte Carlo。在您的情况下，与“1.0”更接近的“x”点对积分值的贡献大于那些更接近“0.0”的点。 - Hristo Iliev

Answer 1

请将以下代码替换：

count = 0
k = 1
while(k <= n){
if(ypoints[k]<=h(xpoints[k]))
    count = count+1
k = k+1
}

使用这行代码：

count <- sum(ypoints <= h(xpoints))

- Matthew Lundberg

---

我认为这相当精辟地概括了向量化的威力。 - Ari B. Friedman

Answer 2

如果你真的追求效率，对于这个问题而言，整合会更快（更节省内存），相比之下要快几个数量级。

integrate(h, 0, 1)

# 1.35 with absolute error < 1.5e-14

microbenchmark(integrate(h, 0, 1), estimate_to(3), times=10)

# Unit: microseconds
#                expr        min         lq     median         uq        max neval
#  integrate(h, 0, 1)     14.456     17.769     42.918     54.514     83.125    10
#      estimate_to(3) 151980.781 159830.956 162290.668 167197.742 174881.066    10

- jbaums