標準偏差、母分散についての信頼区間の求め方とJuliaを用いた計算例

はじめに
分散の信頼区間の計算
データ例を用いたJuliaによる計算
まとめ
参考文献
1. 関連

はじめに

メーカーで働くとき、部品の性能について頭を悩ますことは多いと思う。基本的な測定誤差の解析や考察には計測における誤差解析入門がとても読みやすかったので是非購入してほしい。特に今回は分散値はどのくらいのサンプル数があればある程度信頼できるものになるかを計算する。

より具体的に説明すると、部品というのは典型値(typ)だけではなく、最大/最小値(max/min)も考慮してモジュールを作る必要がある。このとき最小/最大にしきい値を標準偏差σで管理することが多い。例えば3σをしきい値とすると99.7%の歩留まりとなるといったことが予想できるのである。

今回の記事では、このσといった値はどのくらいの個数だとどのくらいのばらつきになるかを計算したので、その計算過程をjuliaで示しながら解説したいと思う。

分散の信頼区間の計算

サンプルの分布が正規分布に従うときに母分散については、下記式はサンプル数\(n\)としたとき\(n-1\)のχ²分布にしたがうことが知られています。
\begin{equation}
\chi^2(n-1) = \frac{(n-1)s^2}{\sigma^2}
\end{equation}
ここで\(s\)は不変分散です。この解説はネットの記事では多くありますが、簡単にイメージをつかむのにはこちら(母分散の意味と区間推定・検定の方法)の記事などがわかりやすいと思います。
これをσについて解いていくと
\begin{equation}
\sigma^2 = \frac{(n-1)s^2}{\chi^2(n-1)}
\end{equation}
ここで\(a\)%信頼区間を考えると以下の範囲に分散は収まる。(例えばa=95のときは95%の確立でこの範囲に分散は収まりますよといったことを表します。)
\begin{equation}
\frac{(n-1)s^2}{\chi^2_{\alpha/2}(n-1)} \leq \sigma^2 \leq \frac{(n-1)s^2}{\chi^2_{1-\alpha/2}(n-1)}
\end{equation}
この式をちゃんと説明すると
左側と右側の値: これが信頼区間の「下限」と「上限」。
\( (n−1)s^2\) : これは手元のサンプルデータから計算できる値。
\(\chi_{α/2}^2(n−1)\) と \(χ_{1-α/2}^2(n−1)\): この値は、サンプルサイズ\(n\)と、どれくらいの確信度\((1-\alpha )\)にしたいかによって決まり、 \(\chi_{α/2}^2(n−1)\)は、χ²分布のグラフで、右側の面積（確率）が \(2\alpha\) になる点です。同様に\(χ_{1-α/2}^2(n−1)\)は、χ²分布のグラフで、左側の面積（確率）が\(2\alpha\)になる点です（右側の面積が \(1-2\alpha\) になる点と同じです）。

データ例を用いたJuliaによる計算

using Random, Distributions, PyPlot

Random.seed!(123)
d = Normal(1000,5)
n = 10000
X = rand(d, n)

# xがちゃんと正規分布をしているかを確認
figure()
x_bins = range(980, 1020, step=1)
hist(X, bins = x_bins)
xlabel("x")
ylabel("frequency")
grid()
display(gcf())

# χ²検定を行う、信頼区間は95%信頼区間
n = size(X)[1]
x⁻ = mean(X)
s² = 1/(n-1)*sum([(xᵢ-x⁻)^2 for xᵢ in X])
γ = 0.95
χₘᵢₙ = quantile(Chisq(n-1) , (1+γ)/2)
χₘₐₓ = quantile(Chisq(n-1) , (1-γ)/2)
σₘᵢₙ = √((n-1)s²/(χₘᵢₙ))
σₘₐₓ = √((n-1)s²/(χₘₐₓ))

println(σₘᵢₙ, "<σ<",σₘₐₓ)

ひとまず、サンプル数n=10000で正規分布と母分散の信頼区間がどのようになるかを計算してみた。この時平均を1000、標準偏差(~√分散)=5の正規分布を仮定している(d = Normal(1000,5))。
このときrand関数で正規分布に従うようにランダムに実験値を出力しており、このときのヒストグラムが以下のようになる。

上記のグラフをみると、きれいに正規分布が再現できていることがわかる。グラフの作成についてはこちらの記事も参考にしてください。
また、このときの標準偏差の95%信頼区間も4.95<σ<5.04となっており、真の値σ=5に近いことがわかる。

ここで、特に自分としては気になったのがサンプル数によって標準偏差の95%信頼区間がどのように変化していくかである。

using Random, Distributions, PyPlot
using PyCall
PyCall.PyDict(matplotlib["rcParams"])["font.family"] = "Meiryo";

Random.seed!(123)
d = Normal(1000,5)
σₘᵢₙ_list = []
σₘₐₓ_list = []

for n in range(2,10000)
    X = rand(d, n)
    # χ²検定を行う、信頼区間は95%信頼区間
    n = size(X)[1]
    x⁻ = mean(X)
    s² = 1/(n-1)*sum([(xᵢ-x⁻)^2 for xᵢ in X])
    γ = 0.95
    χₘᵢₙ = quantile(Chisq(n-1) , (1+γ)/2)
    χₘₐₓ = quantile(Chisq(n-1) , (1-γ)/2)
    σₘᵢₙ = √((n-1)s²/(χₘᵢₙ))
    σₘₐₓ = √((n-1)s²/(χₘₐₓ))

    push!(σₘᵢₙ_list,σₘᵢₙ)
    push!(σₘₐₓ_list,σₘₐₓ)

end

N =  range(2,10000)

function moving_average(x, window_size)
    return [mean(x[max(1, i - window_size + 1):i]) for i in 1:length(x)]
end

# 平滑化：ウィンドウサイズ100くらいで
window_size = 100
σₘᵢₙ_smooth = moving_average(σₘᵢₙ_list, window_size)
σₘₐₓ_smooth = moving_average(σₘₐₓ_list, window_size)

# いろいろと可視化して確認
figure()
plot(N, σₘᵢₙ_list, label="σmin")
plot(N, σₘₐₓ_list, label="σmax")
#ylim(2,5)
legend()
grid(true)
ylabel("分散(σ)")
xlabel("サンプル数")
display(gcf())

figure()
plot(N, σₘᵢₙ_list, label="σmin")
plot(N, σₘₐₓ_list, label="σmax")
ylim(4,6)
legend()
grid(true)
ylabel("分散(σ)")
xlabel("サンプル数")
display(gcf())

figure()
plot(N, σₘᵢₙ_smooth, label="σmin (smoothed)")
plot(N, σₘₐₓ_smooth, label="σmax (smoothed)")
#ylim(2, 3.5)
legend()
grid(true)
xscale("log")
yscale("log")
ylabel("分散(σ)")
xlabel("サンプル数")
xlim(1)
display(gcf())

figure()
plot(N, σₘᵢₙ_smooth, label="σmin (smoothed)")
plot(N, σₘₐₓ_smooth, label="σmax (smoothed)")
ylim(4,6)
legend()
grid(true)
ylabel("分散(σ)")
xlabel("サンプル数")
display(gcf())

この結果が以下のようになる。

これらのグラフからいくつかわかることがわかる。(追記：プログラム中及び図中の分散という表記は間違いで正しくは標準偏差です)
１枚目２枚目からは観測値がランダムなとき、その観測値によってはかなりランダムに信頼区間の幅は変わること。３枚目４枚目はランダムさが大きいところに平滑化処理を加えたものとなっており、σが有効数字一桁の範囲でほしい場合(σ~5)はサンプル数n=100となっており、それ以上に有効数字が必要な時はサンプル数が一気に増えていくことがわかる。