R ビュフォンの針の問題について(確認) R > buffon<-function(n){ L<-1 P<-2 count<-0 for(i in c(1:n)){ x<-runif(1,0,P/2) theta<-runif(1,0,pi/2) if(x<P/2){ if((sin(theta)*L/2)>x) count<-count+1 } else{ if((sin(theta)*L/2)<x) count<-count+1 } } p<-2*L*n/(P*count) p } > buffon(10) [1] 5 > buffon(100) [1] 3.846154 > buffon(1000) [1] 3.215434 > buffon(10000) [1] 3.072197 > buffon(100000) [1] 3.138535 > buffon(1000000) [1] 3.143478 > buffon(10000000) [1] 3.142044 R

`x<-runif(1,0,P/2)` (`x: [0, P/2]`)`if(x

R vector


buffon <- function(n) {
L <- 1; P <- 2
x <- runif(n, 0, P/2)
proj <- sin(runif(n, 0, pi/2))*L/2
2*L*n/(P*sum(proj > x))
}

for(n in 10^c(1:8)) {
cat(sprintf("n = %9d, %f\
", n, buffon(n)))
}

n = 10, 10.000000
n = 100, 2.857143
n = 1000, 3.115265
n = 10000, 3.138732
n = 100000, 3.124219
n = 1000000, 3.140556
n = 10000000, 3.139620
n = 100000000, 3.142656