Monte Carlo pi方法

我嘗試計算R中的Monte Carlo pi函數。我在代碼中遇到了一些問題。現在我寫這篇文章的代碼：Monte Carlo 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

這不是工作，但我不知道如何解決它。

我也想寫一個代碼來繪製這個（用圓形在圓形和圓點內）。

來源

2013-03-03 Phantom

什麼是你想實現什麼？ Monte Carlo pi是什麼？什麼不工作？ – 2013-03-03 13:38:30

在這種情況下，這是用圓形內部正方形計算pi值的方法。我得到這個警告信息： In if（（x^2 + y^2）<= 1）{：條件長度大於1並且只有第一個元素會被使用 – Phantom 2013-03-03 14:22:49

這裏是一個矢量版本（也有出錯了數學）

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")

，但我'd建議你使用較小的N值，也許是10000.

來源

2013-03-03 13:54:23 flodel

+1！真乾淨！ – agstudy 2013-03-03 13:58:21

（+1）請注意，如果太多地增加N，則容易超出RAM。您可以將其包裝在一個循環中（例如，使用sapply）並計算許多估計的平均值，但對於可以實現的精度，該方法非常有限。 – Roland 2013-03-03 14:06:37

順便說一下，它應該是'is.inside < - （x^2 + y^2）<= R^2'，整個事物在'R'中是不變的，所以你可以用1代替。 – Roland 2013-03-03 14:09:17

這是一個有趣的遊戲 - 它有很多版本在網上流動。這裏有一個我從指定來源入侵（他的代碼有點天真）。

從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) 
}

錯字 'n面' 到 '裏面'

來源

2013-03-03 14:26:31

Monte Carlo pi方法

回答

相關問題