Julia Advent Calendar 2023の2日の記事です。この記事はJulia Advent Calendar 2023用に書かれました。

Juliaで確率分布

Juliaで確率分布を扱うにはDistributions.jlを使います。確率分布でよく使われる正規分布は以下のように導入されます。

using Distributions
dist=Normal()

Distributions.Normal{Float64}(μ=0.0, σ=1.0)

Normal型のオブジェクトをdistとして生成します。引数は何も指定していないので、平均1、標準偏差0の正規分布が作られます。確率密度関数のグラフとして表示するにはStatsPlots.jlを使うと簡単です。

using StatsPlots
plot(dist, yrotation=90, label="Normal")
savefig(joinpath(@OUTPUT, "Normal.svg"))
nothing

plot関数にNormal型のオブジェクトを渡すだけで良いです。 savefigでsvgとして保存していますが、これはこのサイトがFranklin.jlで作られており、表示のためにそうしているだけです。

正規分布

正規分布が統計学でよく出てくるのは中心極限定理の所為です。母集団の分布がなんであれ(例外はある)無作為に繰り返し測定すればその測定値の平均値は、母集団の平均値μ(母平均)を中心として分散 $σ^2$ (母分散)をサンプル数nで割った値 $σ^2/n$ の分散を持つ正規分布に近づいていくということです。母集団の分布なんてわからないことが普通なので、独立同分布(iid)の条件を前提におけるのであれば中心極限定理は有用です。 n=1であれば平均μ、分散 $σ^2$ の正規分布となります。

母集団が一様分布の場合

母集団が一様分布の場合を考える。

uniform = Uniform(2,5)
plot(uniform, yrotation=90, label="uniform", ylim=(0,0.5))
savefig(joinpath(@OUTPUT, "uniform.svg"))
nothing

mean(uniform), var(uniform)

(3.5, 0.75)

この母集団の平均は3.5であり、分散は0.75で計算される。この母集団から10個サンプリングする作業を100回繰り返す。

using Random
Random.seed!(1)
samples = [rand(uniform,10) for _ in 1:100]
samples[1:5]

5-element Vector{Vector{Float64}}:
 [2.1475154664443634, 2.3572364492225213, 3.1798130696758418, 2.072282931573583, 4.075571862602665, 4.302554162262173, 2.261759146738227, 4.56715305232872, 4.407682129770471, 3.984276055054304]
 [4.832954761529081, 3.383150168924247, 4.491002993257353, 3.7193952876134206, 2.5298752634385346, 2.344805892886096, 4.3592004955513275, 4.677792704745261, 2.6217587682890557, 2.763417385500673]
 [3.2946421747434216, 3.924646197840194, 3.9031960581667002, 2.9255759467206275, 2.20593736441545, 4.791120961099627, 2.4762826673558416, 4.469367748926269, 2.50115872786336, 4.159162077621555]
 [3.730080597990928, 4.227565336128935, 4.57056890336396, 2.2754446807921274, 2.673783746569811, 4.517452593856156, 2.5058147778099986, 3.6950118846683644, 2.5280515122991853, 3.5921746168734296]
 [3.3757604756601043, 3.925119448843251, 3.9287469451790886, 4.8638801916153, 4.332841322783667, 4.165309850877898, 4.725422833458447, 3.629596002409319, 2.7673697219291262, 3.8494054415521295]

ログは5つまでを表示している。それぞれ平均を計算して、そのヒストグラムを表示する。

sample_means = mean.(samples)
histogram(sample_means, yrotation=90, label="sample_means", bins=20)
savefig(joinpath(@OUTPUT, "histogram.svg"))
nothing

100個程度では正規分布しているようには見えない。平均と分散は

mean(sample_means), var(sample_means)

(3.5315617530433223, 0.08427463601015561)

と計算される。母集団の母平均と母分散をn=10で割った値は

mean(uniform), var(uniform) / 10

(3.5, 0.075)

なので近くはある。

100回の繰り返しは少ないので、5000回繰り返す。

samples = [rand(uniform,10) for _ in 1:5000]
sample_means = mean.(samples)
histogram(sample_means, yrotation=90, label="sample_means", bins=20)
savefig(joinpath(@OUTPUT, "histogram2.svg"))
nothing

正規分布に近づいた。平均と分散は

mean(sample_means), var(sample_means)

(3.4956379305596332, 0.07402950947217028)

であり、こちらも理論上の値とほぼ同等である。