やること
統計学で出てくる95%信頼区間とかいう概念が信用できなかった。
PythonのライブラリであるScipyに統計に関するモジュールがあるので、データから母平均の信頼区間を求めることでPythonで統計を扱う練習をしつつ、信頼区間の真偽を暴く。
環境
Python3.8.5 Jupyter Notebookを用いて行った。
Numpy, Scipyのバージョンはそれぞれ1.19.2, 1.5.2であった。
方法
1. データの作成
numpy.random.normalモジュールを用いてデータの作成を行った。このモジュールは正規分布に従う乱数を返すモジュールであり、平均と標準偏差を引数から指定することができる。
今回はこちらのサイトを参考に、17歳男子の身長の統計を元に擬似的にデータを作成した。
身長の分布は一般に正規分布に従うと考えられている。サイトによると17歳男子の平均身長は170.6931 cm、標準偏差は5.802 cm^2ということがわかったので、モジュールを用いてこれに従う乱数を10人分生成し、擬似的に17歳男子の身長のデータ(仮)を作成することとする。
また、numpy.random.normalモジュールは第３引数にサイズを指定することで、複数のデータを一括に生成することができる。データの作成にあたって便利なのでこれを利用した。
2. t分布を用いてデータの信頼区間を求める
1で作成したデータは一旦天から降ってきたものとし、母平均と母分散は未知であると考える。
データ数が$n$のデータの平均$\mu'$と不偏分散$s^2$を求める。不偏分散とは、データの分散の和を$n-1$で割った値のことである。
自由度$n-1$のt分布における上側2.5%点の値を$t_{n-1}(2.5\%)$と表記する。この値はscipy内のstats.tモジュールを用いることで受け取ることができる。Scipyの統計モジュールstatsで統計分布を使いこなす | 化学の新しいカタチを参考にした。
これらの値から、母平均$\mu$の95%信頼区間は以下のように表せる。
$$\mu'-t_{n-1}(2.5\%) \sqrt{\frac{s^2}{n}} \leq \mu \leq \mu'+t_{n-1}(2.5\%) \sqrt{\frac{s^2}{n}}$$
3. 得られた信頼区間が正しいか確かめる
$\alpha\%$信頼区間とは、同じ方法で繰り返し信頼区間を求めたときに$\alpha\%$の確率で区間内に母平均が含まれるような区間のことである。つまり2の方法で$95\%$信頼区間を求めたとき、100回のうち95回程度は区間内に母平均が含まれているはずである。よって1,2の操作を複数回行い、これが事実かどうか確かめる。
コード
#1. データの作成
data_list = np.random.normal(ave, sd, size)
#データの分散、平均など必要なデータを求める
data_ave = sum(data_list)/size
data_var = sum((data_list-data_ave)**2)/(size-1)
data_t = stats.t.ppf(0.975, size-1)
#2. t分布を用いてデータの信頼区間を求める
#データの平均からの誤差を算出し、信頼区間の上限値、下限値を求める
error = data_t*np.sqrt(data_var/size)
bottom, up = data_ave-error, data_ave+error
1,2の処理を上のコードにまとめた。母平均、母データの標準偏差、データ数となる定数をそれぞれave, sd, sizeとして定義し、計算後、変数bottom, upにそれぞれ信頼区間の下限値、上限値が入る。
以上の処理を関数practice_statsとしてまとめ、関数を1000回実行し、求めたそれぞれの信頼区間の中に母平均が入った回数をカウントするコードを以下にまとめた。
#母平均、標準偏差、データ数を定義
ave, sd, size = 170.6931, 5.802, 10
#求めた信頼区間内に母平均が含まれる回数をカウントし、出力する
count = 0
for i in range(1000):
bottom, up = practice_stat(ave, sd, size)
if bottom <= ave <= up: count += 1
print(count)
結果
上のコードを複数回実行した結果、出力は942、952、945、958のようになった。
結果から1000回のうち950回程度、信頼区間内に母平均が含まれていることがわかる。よって95%信頼区間は文字通り、95%程度信頼できることが事実として理解できた。
より具体的な結果
実際にpractice_stats関数から得られた信頼区間は以下のようになった。
170.5841 <= ave <= 178.3930
166.4153 <= ave <= 175.6009
166.8675 <= ave <= 175.9985
165.9654 <= ave <= 174.9458
165.7519 <= ave <= 175.9660
164.6087 <= ave <= 174.2426
169.8855 <= ave <= 179.0061
164.0898 <= ave <= 173.8109
166.7745 <= ave <= 175.0380
167.4358 <= ave <= 176.5277
(母平均: 170.6931)
10回分の結果をランダムに生成した所、全ての信頼区間内に母平均が含まれていた。95%信頼区間は当然5%の確率で区間内に母平均が含まれないので、何回か実行すると誤った結果が現れることもあった。
数値を変えて実験してみた結果を以下に示す。
(元の母平均: 170.6931、標準偏差: 5.802、データ数: 10)
①標準偏差を5倍にしてみた場合
151.2301 <= ave <= 201.4208
163.1916 <= ave <= 201.3858
164.9185 <= ave <= 196.7287
163.1271 <= ave <= 209.7147
144.3538 <= ave <= 199.3029
155.4677 <= ave <= 197.6442
163.6688 <= ave <= 193.7708
161.7466 <= ave <= 199.4978
134.7038 <= ave <= 174.7072
157.9682 <= ave <= 198.6466
count = 947/1000 (rate = 94.7%)
明らかに信頼区間の幅が広くなっていることがわかる。標準偏差はデータの不偏分散の値に影響するため当然である。
②データ数を1000にしてみた場合
170.2620 <= ave <= 170.9729
170.3597 <= ave <= 171.0572
170.7633 <= ave <= 171.4656
169.9968 <= ave <= 170.7285
170.4933 <= ave <= 171.2160
170.5394 <= ave <= 171.2390
170.0713 <= ave <= 170.7757
170.3891 <= ave <= 171.1146
170.6276 <= ave <= 171.3334
170.4956 <= ave <= 171.1965
count = 964/1000 (rate = 96.4%)
精度は変わらず、信頼区間の幅が非常に狭くなっていることがわかる。3番目の区間の下限値が惜しい。
③99.9%信頼区間にした場合
167.4883 <= ave <= 182.4917
161.2694 <= ave <= 184.6227
158.3516 <= ave <= 179.8887
165.5846 <= ave <= 177.8863
160.9575 <= ave <= 180.8796
162.7626 <= ave <= 181.6215
161.7645 <= ave <= 179.9075
162.6220 <= ave <= 180.0634
162.1130 <= ave <= 184.4436
166.4731 <= ave <= 181.3969
count = 1000/1000 (rate = 100.0%)
信頼区間の幅を犠牲にした代わりに非常に信頼できる区間になった。
scipy.stats.t.interval
実は信頼区間を求めるためにわざわざ式を実装する必要はなく、scipyのstats.t.intervalで信頼度%、自由度、標本平均、標準誤差を指定することで下限値と上限値を受け取る事が可能。
今回の記事では学習の整理が目的だったため、利用しなかった。
結論
Pythonで統計学勉強のモチベを上げることが可能。
参考文献
Scipyの統計モジュールstatsで統計分布を使いこなす | 化学の新しいカタチ
【例】身長の分布は本当に正規分布に従うのか！？ | AVILEN AI Trend
scipy.stats.t — SciPy v1.8.0 Manual
統計学の時間 | 統計WEB