julia で t 分布に従う 95% 信頼区間を求めてみよう

お題のとおり、julia で t 分布に従う平均の 95% 信頼区間を求めてみます。

Distributions.jl を使います。少々、矛盾しますが、100 個の正規分布の乱数を観測値 (ベクトル) とします。 

julia> using Distributions

julia> x = randn(100)
100-element Array{Float64,1}:
 -0.307257
 -0.074242
  0.485464
  1.47269
 -0.0526559
 -1.533
 -0.123655
  0.3673
 -1.62462
 -1.09255
  ⋮
 -0.313221
 -0.655702
 -0.933212
  0.553407
 -1.17567
 -0.305231
  1.8448
  0.0459129
  1.88275
 -0.764762

julia> using PyCall, PyPlot

julia> @pyimport seaborn as sns

julia> sns.distplot(x, kde= false)
PyObject <matplotlib.axes._subplots.AxesSubplot object at 0x130be8940>

x = randn(100) は、x2 = rand(Normal(),100) としても良いです。 

julia> x2 = rand(Normal(),100);

julia> sns.distplot(x2, kde= false)
PyObject <matplotlib.axes._subplots.AxesSubplot object at 0x130be8940>

薄い茶色が後の x2 の方です。

自由度 99 の T 分布 TDist(99) における 2.5% と、97.5% の確率密度関数の x 軸の値は、 quantile() を使って求めることができます。 

julia> quantile.(TDist(99), [0.025, 0.975])
2-element Array{Float64,1}:
 -1.98422
  1.98422

[0.025, 0.975] という2つの値を同時に計算したいので、quantile.() と "." を入れてブロードキャストというのにしています。

もし、これが、正規分布 Normal() だと、1.96 という有名な数字が出てきます。 

julia> quantile.(Normal(), [0.025, 0.975])
2-element Array{Float64,1}:
 -1.95996
  1.95996

やはり、t 分布だと少し広くなりますね。 

julia> mean(x) # 平均
-0.020108918166004477

julia> mean(x) + quantile.(TDist(99), [0.025, 0.975]).* std(x)/(sqrt(100-1))
2-element Array{Float64,1}:
 -0.222323
  0.182105

t 分布の 95% 信頼区間は、 -0.222323 と 0.182105 の間ということになります。

では、観測点を 1,000 個にしてみましょう。 

julia> x3 = randn(1000);

julia> mean(x3)
-0.031912429047342365

julia> mean(x3) + quantile.(TDist(999), [0.025, 0.975]).* std(x3)/(sqrt(1000-1))
2-element Array{Float64,1}:
 -0.0938459
  0.030021

julia> sns.distplot(x3, kde=false, norm_hist=true )
PyObject <matplotlib.axes._subplots.AxesSubplot object at 0x13b4ba198>

1,000個でも滑かとは、言い難いですね。自分達のしている観測って、ちょっと、不安になりますね。

B! LINE