3

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

これは機能していませんが、修正方法がわかりません。

これをプロットするコードも書きたいと思います（四角の中の円とドット）。

4

2 に答える 2

9

これがベクトル化されたバージョンです（そしてあなたの数学にも何か問題がありました）

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を使用することをお勧めします。

于 2013-03-03T13:54:23.870 に答える
0

これは楽しいゲームです-そしてそれの多くのバージョンがウェブの周りに浮かんでいます。これが私が名前の付いたソースからハッキングしたものです（彼のコードはややナイーブでした）。

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'

于 2013-03-03T14:26:31.957 に答える