71772


Rを用いてビュフォンの針の問題をプログラミングしたのですが、実行したところ正しいのかどうかわかりませんでした。（おそらく正しいとは思うものの、確認させていただきます。）

とりあえず、私は以下のようにプログラムをしました。

[code]

    > 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
    }
    
[/code]

実行結果

[code]

    > 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
    
[/code]

ビュフォンの針の問題で求めるのは、円周率の近似値なので、正しいとは思うのですが、Rを始めてからまだ一か月足らずの未熟者で、不安なので一応確認させていただきます。  
よろしくお願いします。


CC BY-SA 4.0
2020-11-06T07:25:25.153
71768
2020-11-06T11:51:34.883
2020-11-06T08:21:31.560
3060
42623
question
r
R ビュフォンの針の問題について(確認)


`x<-runif(1,0,P/2)` としていますので(`x: [0, P/2]`)、`if(x<P/2)` の分岐は必要ないかと思います。  
ところで、R では vector での計算が可能ですので、以下の様にも書くことができます。

[code]

    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", 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
    
[/code]


CC BY-SA 4.0
2020-11-06T10:15:36.333
71772
2020-11-06T11:51:34.883
2020-11-06T11:51:34.883
71768
answer