統計学の備忘録的な学習メモ（1） - 工場裏のアーカイブス

現在、故あってPython（Scipyなど）を使って統計学を基礎から勉強しております。その中で特に備忘録として残しておきたいことをメモして行こうと思います。個人的な備忘録の意味合いが強く、他記事と書き方が異なっており、説明なども少なめです。今後は予告なく内容を追加・修正するかもしれません。

※環境はJupyter Notebookの想定です。
　

インポートするライブラリ

import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
from scipy.special import gamma

from matplotlib import pyplot as plt
import seaborn as sns
sns.set()

さまざまな統計量の表示

Scipyの関数を用いて、(1, 2, 3, 4, 5）という単純なサンプルの様々な統計量を表示する。

sample = np.array([1, 2, 3, 4, 5])

print(f'サンプルサイズ（N） = {len(sample)}') 
print(f"合計 = {sp.sum(sample)}") 
print(f"平均 = {sp.mean(sample)}") 
print(f"最大値 = {sp.amax(sample)}")
print(f"最小値 = {sp.amin(sample)}") 
print(f"中央値 = {sp.median(sample)}") 
print(f"第1四分位数 = {stats.scoreatpercentile(sample, 25)}") 
print(f"第3四分位数 = {stats.scoreatpercentile(sample, 75)}") 
print(f"標本分散（分母が N） = {sp.var(sample, ddof=0)}") 
print(f"不偏分散（分母が N-1） = {sp.var(sample, ddof=1)}") 
print(f"標本標準偏差 = {sp.std(sample, ddof=0):.4f}") 
print(f"不偏分散の平方根 = {sp.std(sample, ddof=1):.4f}")

サンプルサイズ（N） = 5
合計 = 15
平均 = 3.0
最大値 = 5
最小値 = 1
中央値 = 3.0
第1四分位数 = 2.0
第3四分位数 = 4.0
標本分散（分母が N） = 2.0
不偏分散（分母が N-1） = 2.5
標本標準偏差 = 1.4142
不偏分散の平方根 = 1.5811

分散の関数を自前で実装して、比較してみる

#自前の標本分散関数
def my_varp(x):
    mu = sp.mean(x)
    v = sp.sum((x - mu) ** 2) / len(x)
    return v

#自前の不偏分散関数
def my_var(x):
    mu = sp.mean(x)
    v = sp.sum((x - mu) ** 2) / (len(x) - 1)
    return v

print(f"標本分散（自前） = {my_varp(sample)}")
print(f"不偏分散（自前） = {my_var(sample)}")
print(f"標本標準偏差（自前）= {np.sqrt(my_varp(sample)):.4f}")
print(f"不偏分散の平方根（自前）= {np.sqrt(my_var(sample)):.4f}")

自前の関数による計算結果は、確かにScipyと一致している。

標本分散（自前） = 2.0
不偏分散（自前） = 2.5
標本標準偏差（自前）= 1.4142
不偏分散の平方根（自前）= 1.5811

Pandasでも基礎的な統計量の表示は出来る。describe関数を用いると要約統計量を表示してくれる（表示される標準偏差（std）は、実際は不偏分散の平方根である模様）

sr = pd.Series([1, 2, 3, 4, 5])

print(f"合計（pd版） = {sr.sum()}")
print(f"平均（pd版）= {sr.mean()}")
print(f"標本分散（pd版） = {sr.var(ddof=0)}")
print(f"不偏分散（pd版） = {sr.var(ddof=1)}")
print(f"標本標準偏差（pd版） = {sr.std(ddof=1):.4}")
print(f"標準偏差平方根（pd版） = {sr.std(ddof=1):.4}")
print("-----")
sr.describe()

合計（pd版） = 15
平均（pd版）= 3.0
標本分散（pd版） = 2.0
不偏分散（pd版） = 2.5
標本標準偏差（pd版） = 1.4142
不偏分散の平方根（pd版） = 1.5811
-----
count 5.000000
mean 3.000000
std 1.581139
min 1.000000
25% 2.000000
50% 3.000000
75% 4.000000
max 5.000000
dtype: float64

不偏分散の平方根に関する混乱

母集団から抽出した、サンプルサイズ $n$ の標本 $(x_1, x_2, …, x_n)$ がある場合

標本に対する統計量

標本平均	$\bar{x} = \frac{ \sum_{i=1}^n x_i }{n}$
標本分散	$s^2 = \frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n}$
標本標準偏差	$s = \sqrt{\frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n}}$

母集団の母数を標本から推定する統計量

平均	※標本平均がそのまま、母平均の不偏推定量となる
不偏分散	$\hat{\sigma}^2 = \frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n-1}$
不偏分散の平方根	$\hat{\sigma} = \sqrt{\frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n-1}}$

この中でも不偏分散の平方根については、最初に手にした統計学の教科書では「不偏標準偏差」と呼ばれていました。しかし別の教科書では、これは母標準偏差に対する不偏推定量ではないので、不偏標準偏差という言い方はしないと書かれており、非常に混乱しました。

ネットでも色々調べたところ、この辺りの統計量の呼び方や、どんな記号を割り当てるかは非常に混乱しているようです。自分でもシミュレーション（後で掲載）をやってみて、これが母標準偏差に対する不偏推定量ではないことが腑に落ちましたので、不偏分散の平方根とそのまま呼ぼうと思います。非常に参考となった外部サイトをいくつか示します。

不偏標準偏差とは?：統計検定を理解せずに使っている人のために
 不偏分散とは？n-1で割る理由を簡単にわかりやすくエクセルの関数を解説！｜いちばんやさしい、医療統計
 不偏分散の平方根は標準偏差の不偏推定量か | ブログ | 統計WEB
　

サンプルサイズと各統計量のシミュレーション

サンプルサイズを増やして行ったときに、各統計量がどのように変化するかをシミュレーションする

#必要なライブラリはインポート済みとする

#乱数のシードを固定
np.random.seed(1)

#平均50, 標準偏差10、分散100の正規分布を利用
loc = 50
scale = 10
population = stats.norm(loc = loc, scale = scale)

#シミュレーション対象とするサンプルサイズ
sample_size = np.array([10, 100, 1000, 10000])

sample_mean_array = np.zeros(len(sample_size)) #平均を格納
sample_varp_array = np.zeros(len(sample_mean_array)) #標本分散を格納
sample_var_array = np.zeros(len(sample_mean_array)) #不偏分散を格納
sample_stdp_array = np.zeros(len(sample_mean_array)) #標本標準偏差を格納
sample_std_array = np.zeros(len(sample_mean_array)) #不偏分散の平方根を格納

#試行回数は1回として、各サンプルサイズにおける統計量をシミュレートする
for i in range(0, len(sample_size)):
    sample = population.rvs(size = sample_size[i])
    sample_mean_array[i] = sp.mean(sample)
    sample_varp_array[i] = sp.var(sample, ddof=0)
    sample_var_array[i] = sp.var(sample, ddof=1)
    sample_stdp_array[i] = sp.std(sample, ddof=0)
    sample_std_array[i] = sp.std(sample, ddof=1)

#各サンプルサイズにおける統計量を表示
for i in range(0, len(sample_size)):
    print(f'サンプルサイズ {sample_size[i]}')
    print(f'・母平均 = {loc}')
    print(f'・母分散 = {scale ** 2}')
    print(f'・母標準偏差 = {scale}')
    print(f'・母標準誤差 = {scale/sp.sqrt(sample_size[i]):.5}')
    print(' ------')
    print(f'・標本平均 = {sample_mean_array[i]:.5}')
    print(f'・標本分散 = {sample_varp_array[i]:.5}')
    print(f'・標本標準偏差 = {sample_stdp_array[i]:.5}')
    print(f'・標本標準誤差 = {sp.sqrt(sample_varp_array[i]/sample_size[i]):.5}')
    print(' ------')
    print(f'・不偏分散  = {sample_var_array[i]:.5}')
    print(f'・不偏分散の平方根  = {sample_std_array[i]:.5}')
    print(f'・不偏標準誤差（不偏分散の平方根/√N） = {sample_std_array[i]/sp.sqrt(sample_size[i]):.5}')
    print('\n')

試行1回のためか、サンプルサイズ100でも不偏推定量と母数のズレが大きい。サンプルサイズ1000となればかなり近い値となる。

サンプルサイズ 10
・母平均 = 50
・母分散 = 100
・母標準偏差 = 10
・母標準誤差 = 3.1623
------
・標本平均 = 49.029
・標本分散 = 141.82
・標本標準偏差 = 11.909
・標本標準誤差 = 3.766
------
・不偏分散 = 157.58
・不偏分散の平方根 = 12.553
・不偏標準誤差（不偏分散の平方根/√N） = 3.9697

サンプルサイズ 100
・母平均 = 50
・母分散 = 100
・母標準偏差 = 10
・母標準誤差 = 1.0
------
・標本平均 = 50.743
・標本分散 = 68.86
・標本標準偏差 = 8.2982
・標本標準誤差 = 0.82982
------
・不偏分散 = 69.556
・不偏分散の平方根 = 8.34
・不偏標準誤差（不偏分散の平方根/√N） = 0.834

サンプルサイズ 1000
・母平均 = 50
・母分散 = 100
・母標準偏差 = 10
・母標準誤差 = 0.31623
------
・標本平均 = 50.274
・標本分散 = 98.569
・標本標準偏差 = 9.9282
・標本標準誤差 = 0.31396
------
・不偏分散 = 98.668
・不偏分散の平方根 = 9.9332
・不偏標準誤差（不偏分散の平方根/√N） = 0.31411

サンプルサイズ 10000
・母平均 = 50
・母分散 = 100
・母標準偏差 = 10
・母標準誤差 = 0.1
------
・標本平均 = 50.075
・標本分散 = 100.2
・標本標準偏差 = 10.01
・標本標準誤差 = 0.1001
------
・不偏分散 = 100.21
・不偏分散の平方根 = 10.011
・不偏標準誤差（不偏分散の平方根/√N） = 0.10011

標本分布のプロット、および不偏分散の平方根と母標準偏差のズレ

母集団からサイズ10の標本抽出を10,000回繰り返すシミュレーションを行い、標本分布をプロットしてみる。更に、不偏分散の平方根は母標準偏差とどの程度ズレが生じるのかを確認して、参考サイトにある係数c4による補正を試してみる。

#平均50, 標準偏差10、分散100の正規分布を利用
loc = 50
scale = 10
population = stats.norm(loc = loc, scale = scale)

sample_size = 10
n_trial = 10000
sample_mean_array = np.zeros(n_trial)
sample_var_array = np.zeros(len(sample_mean_array))
sample_std_array = np.zeros(len(sample_mean_array))

np.random.seed(2)
for i in range(0, n_trial):
    sample = population.rvs(size = sample_size)
    sample_mean_array[i] = sp.mean(sample)
    sample_var_array[i] = sp.var(sample, ddof=1)
    sample_std_array[i] = sp.std(sample, ddof=1)


#標本平均のヒストグラムに、理論的な標本分布を重ね合わせて表示
sns.distplot(sample_mean_array, kde=False, norm_hist = True)
x = np.arange(start = 40, stop = 60, step = 0.1)
plt.plot(x, stats.norm.pdf(x=x, loc=loc, scale=scale/sp.sqrt(sample_size)))
plt.xlabel("Sample mean")
plt.ylabel("Frequency")

#不偏分散の平方根から、不偏標準偏差（母標準偏差の不偏推定量）を求めるための補正係数
#※ただし、母集団が正規分布に従うことが前提
c4 = sp.sqrt(2/(sample_size-1))*gamma(sample_size/2)/gamma((sample_size-1)/2) 

#各統計量を表示
print(f'サンプルサイズ {sample_size}')
print(f'・母平均 = {loc}')
print(f'・母分散 = {scale ** 2}')
print(f'・母標準偏差 = {scale}')
print(f'・母標準誤差 = {scale/sp.sqrt(sample_size):.5}')
print('--------')
print(f'・標本平均の平均 = {sp.mean(sample_mean_array):.5}')
print(f'・不偏分散の平均 = {sp.mean(sample_var_array):.5}')
print(f'・不偏分散平方根の平均 = {sp.mean(sample_std_array):.5}')
print(f'・不偏標準誤差の平均 = {sp.mean(sample_std_array/sp.sqrt(sample_size)):.5}')
print('--------')
print(f'・c4 = {c4:.5}')
print(f'・不偏標準偏差（不偏分散平方根/c4） = {sp.mean(sp.sqrt(sample_var_array))/c4:.5}')

標本平均の平均、および不偏分散の平均はそれぞれ母数と非常に近い値となっている一方、不偏分散の平方根の平均は9.7366となり母標準偏差（10）よりも小さい。手許の教科書によると、このシミュレーション条件では不偏分散の平方根の理論値は約9.7266となるとのことであり、近い値が得られている。不偏分散の平方根は母標準偏差の不偏推定量ではないことが確認出来た。

そしてこれを係数c4で補正することで、値は10.01となり母標準偏差と非常に近い値となることが確認出来た。

サンプルサイズ 10
・母平均 = 50
・母分散 = 100
・母標準偏差 = 10
・母標準誤差 = 3.1623
--------
・標本平均の平均 = 49.967
・不偏分散の平均 = 100.22
・不偏分散平方根の平均 = 9.7366
・不偏標準誤差の平均 = 3.079
--------
・c4 = 0.97266
・不偏標準偏差（不偏分散平方根/c4） = 10.01