蒙特卡罗π方法

Question

我正在尝试用R计算蒙特卡罗pi函数，但代码存在问题。目前我的代码如下：

ploscinaKvadrata  <- 0
ploscinaKroga <- 0
n = 1000
for (i in i:n) {
  x <- runif(1000, min= -1, max= 1)
  y <- runif(1000, min= -1, max= 1)
  if ((x^2 + y^2) <= 1) {
    ploscinaKroga  <- ploscinaKroga + 1
  } else {
    ploscinaKvadrata <- ploscinaKvadrata + 1
  }
    izracunPi = 4* ploscinaKroga/ploscinaKvadrata
}

izracunPi

这个没有起作用，但我不知道如何修复它。

我还想编写一个代码来绘制它（圆形在正方形内部并带有点）。

- Phantom

---

你想要达成什么？蒙特卡罗方法求圆周率是什么？还有什么不起作用？ - Paul Hiemstra

在这种情况下，这是计算正方形内圆的pi值的方法。我得到了以下警告信息： 在 if ((x^2 + y^2) <= 1) { : 中， 该条件的长度大于1，只有第一个元素会被使用。 - Phantom

Answer 1

这是一个矢量化版本（还有一些数学问题需要修正）。

N <- 1000000
R <- 1
x <- runif(N, min= -R, max= R)
y <- runif(N, min= -R, max= R)
is.inside <- (x^2 + y^2) <= R^2
pi.estimate <- 4 * sum(is.inside) / N
pi.estimate
# [1] 3.141472

关于绘制点的部分，你可以按照以下方式进行操作：

plot.new()
plot.window(xlim = 1.1 * R * c(-1, 1), ylim = 1.1 * R * c(-1, 1))
points(x[ is.inside], y[ is.inside], pch = '.', col = "blue")
points(x[!is.inside], y[!is.inside], pch = '.', col = "red")

但我建议您使用更小的N值，可能是10000。

- flodel

---

请注意，如果您将N增加得太多，很容易超出RAM限制。您可以使用循环（例如，使用sapply）将其包装起来，并计算许多估计值的平均值，但是这种方法在精度方面非常有限。 - Roland

顺便提一下，应该是 is.inside <- (x^2 + y^2) <= R^2，整个过程在 R 中是不变的，所以你可以只使用 1。 - Roland

@Roland。谢谢，我已经修复了。我喜欢保留R，因为它有助于理解。就像文档一样。 - flodel

Answer 2

这是一个有趣的游戏 -- 在网络上有许多版本。这是我从指定来源中修改过来的一个（尽管他的代码有些幼稚）。

来源: http://giventhedata.blogspot.com/2012/09/estimating-pi-with-r-via-mcs-dart-very.html

est.pi <- function(n){

# drawing in  [0,1] x [0,1] covers one quarter of square and circle
# draw random numbers for the coordinates of the "dart-hits"
a <- runif(n,0,1)
b <- runif(n,0,1)
# use the pythagorean theorem
c <- sqrt((a^2) + (b^2) )

inside <- sum(c<1)
#outside <- n-inside

pi.est <- inside/n*4

 return(pi.est)
}

将“nside”拼写错误更正为“inside”

- Carl Witthoft