モンテカルロ 法 円 周 率 – グラン フロント イタリアン イル ギオットーネ

5)%% 0. 5 yRect <- rnorm(1000, 0, 0. 5 という風に xRect, yRect ベクトルを指定します。 plot(xRect, yRect) と、プロットすると以下のようになります。 (ここでは可視性重視のため、点の数を1000としています) 正方形っぽくなりました。 3. で述べた、円を追加で描画してみます。 上図のうち、円の中にある点の数をカウントします。 どうやって「円の中にある」ということを判定するか? 答えは、前述の円の関数、 より明らかです。 # 変数、ベクトルの初期化 myCount <- 0 sahen <- c() for(i in 1:length(xRect)){ sahen[i] <- xRect[i]^2 + yRect[i]^2 # 左辺値の算出 if(sahen[i] < 0. 25) myCount <- myCount + 1 # 判定とカウント} これを実行して、myCount の値を4倍して、1000で割ると… (4倍するのは2. より、1000で割るのも同じく2. より) > myCount * 4 / 1000 [1] 3. 128 円周率が求まりました。 た・だ・し! 我々の知っている、3. 14とは大分誤差が出てますね。 それは、点の数(サンプル数)が小さいからです。 ですので、 を、 xRect <- rnorm(10000, 0, 0. 5 yRect <- rnorm(10000, 0, 0. 5 と安直に10倍にしてみましょう。 図にすると ほぼ真っ黒です(色変えれば良い話ですけど)。 まあ、可視化はあくまでイメージのためのものですので、ここではあまり深入りはしません。 肝心の、円周率を再度計算してみます。 > myCount * 4 / length(xRect) [1] 3. 1464 少しは近くなりました。 ただし、Rの円周率(既にあります(笑)) > pi [1] 3. モンテカルロ法 円周率 考察. 141593 と比べ、まだ誤差が大きいです。 同じくサンプル数をまた10倍してみましょう。 (流石にもう図にはしません) xRect <- rnorm(100000, 0, 0. 5 yRect <- rnorm(100000, 0, 0. 5 で、また円周率の計算です。 [1] 3. 14944 おっと…誤差が却って大きくなってしまいました。 乱数の精度(って何だよ)が悪いのか、アルゴリズムがタコ(とは思いたくないですが)なのか…。 こういう時は数をこなしましょう。 それの、平均値を求めます。 コードとしては、 myPaiFunc <- function(){ x <- rnorm(100000, 0, 0.