いつも拝読している「完全無欠で荒唐無稽な夢」というブログがあるのですが、そのコメント欄にこんなことが書いてありました。
こんにちは ^^
またお願いなんですが…^^;
ビュフォンの針：平行線の間隔の半分の長さの針が平行線と交わる確率が1/πというものですが…
平行線を同心円に変換した場合，同じように針が同心円と交わる確率ってのは求められるのでしょうか知らん…?
平行線も針も曲げた場合は同じ確率になるわけだから...針だけ直線のままなら…1/πよりも大きくなりそうとは予測できるのですが...計算の方法がよくわかりません…^^;;
例の，モンテカルロ法で近似値だけでも分かればと…Orz…
お手すきのときにでもよろしければお願い致します～m(_ _)m～v
うーん、これはおもしろそうだと考えてみたのですが、コンピュータ計算で任意の実数を等確率に与える乱数を得ることは無理だと気付きました。ただ、これは有限範囲でも近い数字は出せそうなので、やってみました。
針の長さは 2、同心円の間隔は 4 で計算します。
コードはこんな風でいいのですかねえ…（自信なし） 使った言語は Ruby です。
include Math R = 100000.0 N = 100_0000 def get_num(x, y) (sqrt(x ** 2 + y ** 2) / 4).to_i end counter = 0 N.times do r = rand * R Θ = rand * (2 * PI) φ = rand * PI n1 = get_num(r * cos(Θ) + cos(φ), r * sin(Θ) + sin(φ)) n2 = get_num(r * cos(Θ) - cos(φ), r * sin(Θ) - sin(φ)) counter += 1 unless n1 == n2 end puts counter / N.to_f
このコードでは、同心円の中心から半径 R =100000.0 以内の場所に針が落ちます。試行回数は N = 1000000 回です。
この例でやってみると、確率はだいたい 0.318 くらいになりました。R の値を変えても確率は大きくは変わりませんでした。ただ、R が小さいと確率は多少低くなります（R = 100.0 で 0.312 くらい）。1 / π がだいたい 0.3183... なので、同心円の場合も確率は 1 / π になるのではないでしょうか。
なお、ここでは針の中心の位置を同心円の中心からの距離と回転角であたえていますが、このあたえ方は本質的ではありません。針が同心円と交わるかを調べるときに、それに依存していないからです。だから、完全にランダムであれば、正方形領域に針を落としても同じ答えになる筈です。
※追記
2点で交わる場合の見落としがありました。こうかな。
include Math R = 100000.0 N = 100_0000 def get_num(x, y) (sqrt(x ** 2 + y ** 2) / 4).to_i end counter = 0 N.times do r = rand * R Θ = rand * (2 * PI) φ = rand * PI n1 = get_num(x1 = r * cos(Θ) + cos(φ), r * sin(Θ) + sin(φ)) n2 = get_num(x2 = r * cos(Θ) - cos(φ), r * sin(Θ) - sin(φ)) x3 = ( sin(φ) ** 2 * cos(Θ) - sin(φ) * cos(φ) * sin(Θ)) * r y3 = (-sin(φ) * cos(φ) * cos(Θ) + cos(φ) ** 2 * sin(Θ)) * r n3 = get_num(x3, y3) x1, x2 = x2, x1 if x1 < x2 if n1 != n2 counter += 1 elsif n1 == n2 and n2 - n3 == 1 and x3 < x1 and x2 < x3 counter += 1 end end puts counter / N.to_f
これでも確率はほぼ 0.318、誤差の範囲で変わらないですね。
BigDecimal を使って精度を上げてみる。しかし、時間がかかりすぎて充分な試行数が取れない…。
require 'bigdecimal' require 'bigdecimal/math' include Math include BigMath R = 100000.0 N = 10000 F = 12 def get_num(x, y) (sqrt(x ** 2 + y ** 2, F) / 4).to_i end counter = 0 N.times do r = BigDecimal(rand * R, F) Θ = BigDecimal(rand * (2 * PI), F) φ = BigDecimal(rand * PI, F) n1 = get_num(x1 = r * cos(Θ, F) + cos(φ, F), r * sin(Θ, F) + sin(φ, F)) n2 = get_num(x2 = r * cos(Θ, F) - cos(φ, F), r * sin(Θ, F) - sin(φ, F)) x3 = ( sin(φ, F) ** 2 * cos(Θ, F) - sin(φ, F) * cos(φ, F) * sin(Θ, F)) * r y3 = (-sin(φ, F) * cos(φ, F) * cos(Θ, F) + cos(φ, F) ** 2 * sin(Θ, F)) * r n3 = get_num(x3, y3) x1, x2 = x2, x1 if x1 < x2 if n1 != n2 counter += 1 elsif n1 == n2 and n2 - n3 == 1 and x3 < x1 and x2 < x3 counter += 1 end end puts counter / N.to_f
これも確率は 0.318 くらい。
どうやら R を大きくすると、2点で交わる場合はほとんど無視できるようですね。結局ほぼ 1 / π ということですか。
※再追記
中心に大きな穴ができるのをなくしました。メソッド get_num を書き換えています。
def get_num(x, y) l = sqrt(x ** 2 + y ** 2) - 2 return 0 if l <= 0 (l / 4).to_i + 1 end
こうしても R が大きい場合の確率は誤差の範囲で変わりません。R を小さくすると効果がでます。R = 100.0 の場合の確率は 0.318 となって、補正されていることがわかります。また R = 10.0 とさらに小さくしてみると、確率は 0.323 と今度は大きめになり、2点で交わる効果が出てきて、これも期待どおりになります。
