2017年7月31日月曜日

平均が簡単なようでいて難しい

簡単に計算できるかな、というか解析的に求まるのかな、と思っていた平均が実はそれほど簡単でもなかった、というのがあって、以下のような問題が思っていた以上に難しい。


問題:「一辺が長さ 1 の正方形にランダムに2点を置いた場合の2点間の距離の平均はいくらか?」

それほど厳密な解がほしいわけでもないので、サクッとモンテカルロシミュレーションを Julia で作ってみたら、だいたい 0.5214 という解になった。
この 0.5214 って数字、なんだろう?


ちなみに、作ったプログラムは
sum0 = 0.0
T = 10000*10000
for i=1:T
    x1 = rand()
    x2 = rand()
    y1 = rand()
    y2 = rand()
    dist0 = sqrt((x1-x2)^2 + (y1-y2)^2)
    sum0 = sum0 + dist0
    if rem(i, 10000000) == 0
        println("iter = $i, average = $(sum0/i)")
    end
end

で 40 秒ぐらいで計算してみている。
julia> @time include("randdist.jl")
iter = 10000000, average = 0.521513779364118
iter = 20000000, average = 0.521471850234503
iter = 30000000, average = 0.5214132474790024
iter = 40000000, average = 0.5214231135462779
iter = 50000000, average = 0.521431221770643
iter = 60000000, average = 0.5214123606643466
iter = 70000000, average = 0.5213983220702314
iter = 80000000, average = 0.521385077693363
iter = 90000000, average = 0.5213964364889536
iter = 100000000, average = 0.5213935109322017
 41.304951 seconds (600.00 M allocations: 10.431 GB, 1.90% gc time)

ちなみに、1行で計算すると4倍ぐらい速かった。

julia> @time begin T = 10000*10000; x = rand(T,4); ave = sum(sqrt((x[:,1] - x[:,2]).^2 + (x[:,3]-x[:,4]).^2))/T; println("T = $T, ave = $ave"); end
T = 100000000, ave = 0.5213702192902951
  8.407286 seconds (45 allocations: 10.431 GB, 9.60% gc time)




0 件のコメント:

コメントを投稿