Juliaで標準正規分布

今回は、統計学の基本のキの正規分布関係。先人のエントリーに学んで淡々とやっていく。なるべく、juliaらしさを活かした書き方を追求したい。

参考 mWsoftさんの標準正規分布など　#３年半くらい前のjuliaエントリー

標準正規分布は、教科書的には

①　f(x) = \frac {1} { \sqrt {2 \pi} } \exp(- \frac {x^2} {2})

なのだろう。

今回は、juliaで書きやすいように

②　f(x) = \frac {exp(- \frac {x^2} {2})} { \sqrt {2 \pi} }

と変形して、
さらに、

③　\left\{
\begin{array}{ll}
fo(x) = - \frac   {x^2} {2} \\
f(x) = \frac {exp(fo(x))} { \sqrt {2 \pi} }
\end{array}
\right.

と分離して、途中の数値を確認できるようにしておく。

juliaで③を書くと

.py

# 標準正規分布
fo(x) = -x^2 / 2
f(x) = exp(fo(x)) / √(2π)

と書ける。

juliaの演算順位から

.py

julia> √(2π) < √2π
true

であることに注意しておこう。演算順位が自明でない時には、かっこをつけておくことにする。

さて、fo,fのそれぞれの値を、確認しておこう。

std1.jl

# 標準正規分布の函数fを定義
fo(x) = -x^2 / 2
f(x) = exp(fo(x)) / √(2π)

# プリント関数（副作用あり)
P(g,n) = let 
    println("函数名:" , g)
    result = [g(x) for x in -n:n] # 内包表記
    println(result)
end

P(fo, 4)
P(f, 4)

...オレオレjulia用語で恐縮だが、副作用がない数学的なfunctionを函数、副作用あるFUNCTIONを関数と呼ぶことにする。また、pythonと概ね同じリスト内包表記を用いている。

実行結果:

julia.exe std.jl
函数名:fo
[-8.0, -4.5, -2.0, -0.5, 0.0, -0.5, -2.0, -4.5, -8.0]
函数名:f
[0.00013383, 0.00443185, 0.053991, 0.241971, 0.398942, 0.241971, 0.053991, 0.00443185, 0.00013383]

0のときが最大値を取り、y軸に線対称な標準正規分布が表現されているようだ。
函数名をそのままprintできるのはいいね。

グラフを描く

なぜかREPLでしか動かなかった。
導入したバッケージへのパスを切ることができていないのかな。
描画関係は環境依存だろうし、気が向いたら調査することとしたい。

Pkg.add("Plots")
としてインストール．その後、
Pkg.add("GR")
とする。

.jl

julia> using Plots

julia> gr()
Plots.GRBackend()

julia> plot(f)

描画結果:

windows環境にQt周りの環境が構築されるらしい。

『現代数理統計学の基礎』の世界へ

この調子で、いろいろな確率分布をjuliaで書いていきたい。
羅針盤とするのは、新し目な教科書である『現代数理統計学の基礎』（共立出版，2017）。

書籍中の演習問題と回答を以下に公開してくれいるらしいので、演習問題周りを重視で。
https://sites.google.com/site/ktatsuya77/xian-dai