正規分布の標準偏差と信頼区間の計算について

正規分布の標準偏差と3σの信頼区間について調べたことの記録。

背景

ある値の分布が正規分布に従っている場合、その標準偏差をσとすると、平均値±σの範囲にある確率は約68%、平均値±3σの範囲内にある確率は99.7%だという。

この概念は工業製品の製造現場などでよく使われていて、ある製品が検査工程を経て市場に出回ったときに、どれくらいの確率で不良品となるかを見積もるときなどに使われる。
こういうときは3σや6σの水準が使われることが多い。6σになると100万分の3、4とかの比率になる。

こういったことは学校や会社ではそういうものと教わったし、そういうものとだけ理解しているけど、確かめたことはない。

WikiPediaの標準偏差についての項目(標準偏差 - Wikipedia)にも数式があるし、それを解けばわかるのだろうけど、わかった気になるだけの可能性が99.7%になりそう。そこで、今回は真面目に数式を解くことはせずに、計算機の力を借りて確かめてみることにした。

環境

Python 3.11
以下のコードはJupyter Notebookでの作業・表示に最適化したもの。

n=100の正規分布

まず、任意の整数nで正規分布の乱数を生成してヒストグラムを出力する関数を定義する。
乱数を生成するのに random.randn() 関数を使用した。標準正規分布(ガウス分布)に基づいた乱数を生成してくれるもので、平均値が0、標準偏差が1になるように正規化されている。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

#Jupyter Notebook上でグラフを表示するための儀式
%matplotlib inline


#任意の数(n)で正規分布のヒストグラムを作製するための関数を定義
# n: 分布を作製するデータの数
# summary: n数、信頼区間のサマリをコンソールに表示するかどうかを選択(Boolean)

def norm_dist_plot(n,summary):
    
    x=np.random.randn(n) # x:正規分布の乱数を生成 
    
    x_avg=np.mean(x) # 平均値(定義上は0)
    x_std=np.std(x) # 標準偏差(定義上は1)

    
    # 3σ 信頼区間の比率の計算
    usl_3=x_avg+x_std*3
    lsl_3=x_avg-x_std*3
    cnt_3sig=np.sum((x<usl_3)&(x>lsl_3))
    sig3_rate=cnt_3sig/n
    
    norm_dist_plot=[n,x_avg,x_std,usl_3,lsl_3,sig3_rate]

    
    if summary==True:
        
        print("n = ",n)
        print("x_avg = ",x_avg)
        print("x_std = ",x_std)
        print("usl(3σ)= ",usl_3)
        print("lsl(3σ)= ",lsl_3)
        print("Ratio within the ±3σ = ",str(sig3_rate) )
    
    
    # 正規分布のヒストグラムを作製
    plt.hist(x,bins=50,normed=True)
    plt.title("n="+str(n))
    
    #理想的な正規分布を示す補助線を追加する
    y = mlab.normpdf(np.arange(-3,3,0.1), 0, 1)
    plt.plot(np.arange(-3,3,0.1),y,'--')
    
    plt.axvline(usl_3, ls = "--", color = "navy",linewidth = 1)
    plt.axvline(lsl_3, ls = "--", color = "navy",linewidth = 1)
    
    return norm_dist_plot

Notebookに上記のコマンドを実行したあと、次のセルで以下のコマンドを実行する。

# n=100で正規分布のヒストグラムを作製
# True⇒コンソールにサマリを出力
n_100=norm_dist_plot(100,True)

出力されたのが以下のサマリとヒストグラム。
Trueをいれるとn数、平均値や標準偏差、平均値±3σにある値の比率などをサマリとして出力するようにしている。
ヒストグラムをみるとさすがにn=100は少なくて、あまり正規分布のような分布にはみえない。しかも、平均値±3σにある値の比率は100%だと判定されてしまっている。たしかにこの分布ではそうなっているし、仕方がない。
ちなみに、オレンジ色の線は標準正規分布を示すための補助線。理想的にはヒストグラムの山がこの線と重なるようになる。

n=10000の正規分布

n数を100倍にして、n=10000でもう一度実行した。

n_10000=norm_dist_plot(10000,True)

出力されたのが以下。
ヒストグラムが正規分布っぽい山になったし、それを示すオレンジ色の補助線にも概ね重なっている。平均値±3σにある値の比率は0.9972(97.2%)とのことで、ほぼ教科書通りの結果となった。

オマケ

最後に、オマケとしてn=100, 1000, 10000, 100000としたときの、それぞれのヒストグラムを出力して並べてみる。

fig=plt.figure(figsize=(10,6))

plt.subplot(2,2,1)
n_100=norm_dist_plot(100,False) #サマリ出力は省略

plt.subplot(2,2,2)
n_1000=norm_dist_plot(1000,False)

plt.subplot(2,2,3)
n_10000=norm_dist_plot(10000,False)

plt.subplot(2,2,4)
n_100000=norm_dist_plot(100000,False)

出力したのが以下のヒストグラム。 n数が大きくなるにつれて、ヒストグラムの山が理想的な正規分布(オレンジ色の線)に近づいていっている様子がわかる。

このとき、サマリを出力しないように、定義した関数 norm_dist_plot()の2項目をFalseにしている。 それぞれのサマリの値はn_100〜n_100000に格納されるようにしているので、それらを並べて比較できるようにする。Pandasライブラリを使用した。

import pandas as pd

a=np.vstack((n_100,n_1000,n_10000,n_100000))
df = pd.DataFrame(a, columns=['n数','平均値','標準偏差','3σ上限値','3σ下限値','平均値±3σ比率'])

df

出力されたのが以下のテーブル。n数が100000になると平均値±3σにある値の比率はほぼ99.7%で、ほぼ教科書どおりとなった。

まとめ

数式をといて証明はしなかったけど、計算機に計算させてそれらしい結果が出たので「ああ、やっぱりそうなんだ……」という気持ちになった。