文系プログラマが数式アレルギーを治す為に、中学〜高校レベルの知識でもわかる数式を、プログラムのコードに直す作業を黙々と続ける企画。
ソースコードはJuliaを利用。「PythonやRよりずっとはやい！」ことで有名だけど、言語としての性能が実行速度の決定的な差ではないことも同時に教えられるJuliaを利用。向こうは性能の為ならFortran持ち出してきたりするから仕方ないね。
前回は標準正規分布をやったので、今回は正規分布。
式はWikipediaの正規分布を参考にした。
http://ja.wikipedia.org/wiki/%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83
こんな感じ。
# 自分用メモ
$ f(x) = \frac {1} { \sqrt {2 \pi \sigma^2} } \exp(- \frac {(x - \mu)^2} {2 \sigma^2}) $
標準正規分布の時と比べて、σとμが増えてるけど、大きく違わないのでそれほど面倒ではなさそう。
π（パイ）は円周率。
σ（シグマ）は標準偏差を意味することが多い。
μ（ミュー）は平均を意味することが多い。この場合は算術平均。
f(x)はプログラマ的には、xを引数に取る関数として考えると馴染み深い。
expは指数関数。exponential function。ネイピア数を底と取るやつ。
μ（平均）とσ（標準偏差）は前にやったので、meanとstdで求めることにする。
a = [3.0, 2.5, 1.4, 6.2, 3.3, 2.1, 5.0] μ = mean( a ) #=> 3.3571428571428577 σ = std( a ) #=> 1.6860774427223513
こうやってμやσを気軽に変数名にできる言語って良いよね。あとは数式通りにコードに変換するだけ。
まずは前半、expより前のところだけ書いてみる。
1 / sqrt( 2 * π * σ ^ 2 )
次にexpの中。
- (x - μ) ^ 2 / 2 * σ ^ 2
これを合わせると、こんな風になる。
f = x -> ( 1 / sqrt( 2 * π * σ ^ 2 ) ) * exp( - ( (x - μ) ^ 2 ) / ( 2 * σ ^ 2 ) )
……合ってるよね？ よくこの辺、書き間違えて変な結果を出してしまうことがある。括弧の数はもう少し減らせるはずだけど、無駄に付けてしまっている。
式が正しければ、平均から標準偏差×±1.96の範囲がおよそ95%になるはず。
ということで標準正規分布の時にやったのと同じ感じで、確率を求めてみる。
d = 1e5 y = 0.0 for i = 0:(d * σ * 1.96) y += f(μ + i / d) * (1 / d) end y * 2 #=> 0.9500067984803418
おっけー。予定通り95%。