歯 の 噛み 合わせ 治し 方 割り箸
Pythonでモンテカルロ法を使って円周率の近似解を求めるというのを機会があってやりましたので、概要と実装について少し解説していきます。 モンテカルロ法とは モンテカルロ法とは、乱数を用いてシミュレーションや数値計算を行う方法の一つです。大量の乱数を生成して、条件に当てはめていって近似解を求めていきます。 今回は「円周率の近似解」を求めていきます。モンテカルロ法を理解するのに「円周率の近似解」を求めるやり方を知るのが一番有名だそうです。 計算手順 円周率の近似値を求める計算手順を以下に示します。 1. 「1×1」の正方形内にランダムに点を打っていく (x, y)座標のx, yを、0〜1までの乱数を生成することになります。 2. モンテカルロ法 円周率 求め方. 「生成した点」と「原点」の距離が1以下なら1ポイント、1より大きいなら0ポイントをカウントします。(円の方程式であるx^2+y^2=1を利用して、x^2+y^2 <= 1なら円の内側としてカウントします) 3. 上記の1, 2の操作をN回繰り返します。2で得たポイントをPに加算します。 4.
024\)である。 つまり、円周率の近似値は以下のようにして求めることができる。 N <- 500 count <- sum(x*x + y*y < 1) 4 * count / N ## [1] 3. 24 円周率の計算を複数回行う 上で紹介した、円周率の計算を複数回行ってみよう。以下のプログラムでは一回の計算においてN個の点を用いて円周率を計算し、それを\(K\)回繰り返している。それぞれの試行の結果を に貯めておき、最終的にはその平均値とヒストグラムを表示している。 なお、上記の計算とは異なり、第1象限の1/4円のみを用いている。 K <- 1000 N <- 100000 <- rep(0, times=K) for (k in seq(1, K)) { x <- runif(N, min=0, max=1) y <- runif(N, min=0, max=1) [k] <- 4*(count / N)} cat(sprintf("K=%d N=%d ==> pi=%f\n", K, N, mean())) ## K=1000 N=100000 ==> pi=3. 141609 hist(, breaks=50) rug() 中心極限定理により、結果が正規分布に従っている。 モンテカルロ法を用いた計算例 モンティ・ホール問題 あるクイズゲームの優勝者に提示される最終問題。3つのドアがあり、うち1つの後ろには宝が、残り2つにはゴミが置いてあるとする。優勝者は3つのドアから1つを選択するが、そのドアを開ける前にクイズゲームの司会者が残り2つのドアのうち1つを開け、扉の後ろのゴミを見せてくれる。ここで優勝者は自分がすでに選んだドアか、それとも残っているもう1つのドアを改めて選ぶことができる。 さて、ドアの選択を変更することは宝が得られる確率にどの程度影響があるのだろうか。 N <- 10000 <- floor(runif(N) * 3) + 1 # 宝があるドア (1, 2, or 3) <- floor(runif(N) * 3) + 1 # 最初の選択 (1, 2, or 3) <- floor(runif(N) * 2) # ドアを変えるか (1:yes or 0:no) # ドアを変更して宝が手に入る場合の数を計算 <- (! モンテカルロ法で円周率を求めるのをPythonで実装｜shimakaze_soft｜note. =) & () # ドアを変更せずに宝が手に入る場合の数を計算 <- ( ==) & () # それぞれの確率を求める sum() / sum() ## [1] 0.
01 \varepsilon=0. 01 )以内にしたい場合, 1 − 2 exp ( − π N ⋅ 0. 0 1 2 12) ≥ 0. 9 1-2\exp\left(-\frac{\pi N\cdot 0. 01^2}{12}\right)\geq 0. 9 ならよいので, N ≒ 1. 1 × 1 0 5 N\fallingdotseq 1. モンテカルロ法と円周率の近似計算 | 高校数学の美しい物語. 1\times 10^5 回くらい必要になります。 誤差 %におさえるために10万個も点を打つなんてやってられないですね。 ※Chernoffの不等式については, Chernoff bounds, and some applications が詳しいです。ここでは,上記の文献の Corollary 5 を使いました。 「多分うまくいくけど失敗する可能性もあるよ〜」というアルゴリズムで納得しないといけないのは少し気持ち悪いですが,そのぶん応用範囲が広いです。 ◎ 確率・統計分野の記事一覧
5 y <- rnorm(100000, 0, 0. 5 for(i in 1:length(x)){ sahen[i] <- x[i]^2 + y[i]^2 # 左辺値の算出 return(myCount)} と、ただ関数化しただけに過ぎません。コピペです。 これを、例えば10回やりますと… > for(i in 1:10) print(myPaiFunc() * 4 / 100000) [1] 3. 13628 [1] 3. 15008 [1] 3. 14324 [1] 3. 12944 [1] 3. 14888 [1] 3. 13476 [1] 3. 14156 [1] 3. 14692 [1] 3. 14652 [1] 3. 1384 さて、100回ループさせてベクトルに放り込んで平均値出しますか。 myPaiVec <- c() for(i in 1:100) myPaiVec[i] <- myPaiFunc() * 4 / 100000 mean(myPaiVec) で、結果は… > mean(myPaiVec) [1] 3. 141426 うーん、イマイチですね…。 あ。 アルゴリズムがタコだった(やっぱり…)。 の、 if(sahen[i] < 0. 25) myCount <- myCount + 1 # 判定とカウント ここです。 これだと、円周上の点は弾かれてしまいます。ですので、 if(sahen[i] <= 0. 25) myCount <- myCount + 1 # 判定とカウント と直します。 [1] 3. 141119 また誤差が大きくなってしまった…。 …あんまり関係ありませんでしたね…。 といっても、誤差値 |3. 141593 - 3. 141119| = 0. モンテカルロ法 円周率 c言語. 000474 と、かなり小さい(と思いたい…)ので、まあこんなものとしましょう。 当然ですけど、ここまでに書いたコードは、実行するたび計算結果は異なります。 最後に、今回のコードの最終形を貼り付けておきます。 --ここから-- x <- seq(-0. 5, length=1000) par(new=T); plot(x, yP, xlim=c(-0. 5)) myCount * 4 / length(xRect) if(sahen[i] <= 0. 25) myCount <- myCount + 1 # 判定とカウント} for(i in 1:10) print(myPaiFunc() * 4 / 100000) pi --ここまで-- うわ…きったねえコーディング…。 でもまあ、このコードを延々とCtrl+R 押下で図形の描画とπの計算、両方やってくれます。 各種パラメータは適宜変えて下さい。 以上!