R中的蒙特卡洛方法答案

【问题标题】：Monte Carlo Method in RR中的蒙特卡洛方法
【发布时间】：2018-11-12 20:53:21
【问题描述】：

我正在尝试学习 R。我正在尝试编写一个计算（大约）pi 的程序。 Read About the method

我的代码现在不工作！

f <- 0
s <- 0
range <- 10000
for (i in (1:range)) {
    v <- sample(1:range, 1)/range
    n <- sample(1:range, 1)/range

    if ( sqrt (v*v + n*n) <= 1) {
        f <- f + 1
    } else if ( v <=1 && n <= 1) {
        s <- s+1
    }
}

print ( f/s )

【问题讨论】：

在给定的代码中，f 位于圆圈内，s 由于 else if 语句而位于外部圆圈内。因此，f/s 的结果似乎是 (pi/4)/(1-pi/4) = 3.6597... 如果只用 if 替换 else if，那么结果是 (pi/4) /1，乘以 4 后得到 4*(pi/4) = pi。
...实际上，你也可以删除if( v <= 1 && n <= 1)，因为sample(1:range, 1)/range总是在1/range和1之间。为了符合维基百科的方法，你需要sample(0:range, 1)/range，即使在实践中（range 的大值）也不会太重要
最后，您的代码适合 R 中的第一步，但请注意，R 与 C 或 Java 不同，因为您需要对所有内容进行循环。大多数情况下，您可以一次修改整个向量，这更快，更重要的是，更短（更少的错误）。@mickey 的解决方案表明了这一点。例如，他将一个整数向量与给出逻辑向量的 1 进行比较。 sum 然后只计算该向量中的 TRUEs。随着时间的推移，我建议您通过R Inferno 工作。例如，您会在第 3 圈中烧毁（“无法矢量化”）；-)

标签： r montecarlo

【解决方案1】：

这是您的代码的改进版本

range = 100000
v = runif(range)
n = runif(range)
f = sum(sqrt(v^2 + n^2) <= 1)

print(4 * f / range)

您应该使用runif 而不是sample(...) / range 从制服中获取样本。

s 是不必要的，因为您正在做的是计算次数，f，您的随机点 (v,n) 在该象限的圆圈内，除以尝试抽奖的次数，即在你的情况下只是range。

您需要乘以4，因为f / range 近似于单位圆四分之一的面积。

【讨论】：