Pythonで学ぶあたらしい統計学の教科書
のまとめやコードを実行してみた際のメモです．
試行，標本分布
試行：一回の調査や実験のこと．全く同じ条件下で調査や実験を繰り返すことができる場合は，試行を複数回行うことができる．
試行回数：何度も試行を繰り返すことができる場合，繰り返し回数のことを試行回数と呼ぶ．
標本分布：標本の統計量が従う確率分布のこと．
通常は，調査は1回のみしか行いません（全く同じ調査を，同じ条件で行うことは，ほぼない...と言って良い）．サンプルサイズに関わらず，1回の調査から得られる標本は1つのみです．以下で考える標本分布は「もしも」同じ調査を何度も繰り返したならばどうなるか（例えば，試行回数が3回の場合は，標本平均が３つ得られる）ということを考えたものです．
まずは，数値計算に必要なライブラリをインポートします．
>>> import numpy as np
>>> import pandas as pd
>>> import scipy as sp
>>> from scipy import stats
グラフを描画するライブラリをインポートします．
>>> from matplotlib import pyplot as plt
>>> import seaborn as sns
>>> sns.set()
REPLでは必要ありませんが，Jupyter Notebookでは表示桁数の指定と，グラフを表示させるための指定を行います．
>>> %precision 3
>>> %matplotlib inline
ここでは，母集団は常に平均 4.00, 標準偏差 0.80（分散 0.64）の正規分布とすることにします．これらの条件を設定した population を定義します．
>>> population = stats.norm(loc = 4, scale = 0.8)
標本平均の計算
試行回数を10000回としてみます（標本平均が10000個得られる）．
平均値を格納する配列（長さ 10000のアレイ）を作成し，10000個の標本平均を格納します．
>>> sample_mean_array = np.zeros(10000)
>>> np.random.seed(1)
>>> for i in range(0, 10000):
... sample = population.rvs(size = 10)
... sample_mean_array[i] = sp.mean(sample)
1行目：乱数の種を指定しています．
2行目：for構文を用いて10000回の繰り返しを指定しています．
3行目以降：母集団からの標本抽出シミュレーション．sample_mean_array のi番目に，計算された標本平均を格納しています．標本の抽出は population.rvs 関数で行い（stats.norm.rvs(loc = 4, scale = 0.8, size = 10 としても同じこと），サンプルサイズは10になります．得られた標本に対して sp.mean 関数を用いて平均値を計算しています．
10000個の標本平均は以下のようになります．
>>> sample_mean_array
array([3.92228729, 3.86432929, 4.06953003, ..., 4.13616562, 4.57489661,
4.09896685])
標本平均の平均値と母平均の関係
標本平均の平均値を計算します．シミュレーションで 10000個の標本平均が得られており，10000個の標本平均の平均値を得るという計算になります．
>>> sp.mean(sample_mean_array)
4.004202422791747
母平均は 4 なので，かなり近い値になっています．
標本平均の標準偏差を計算してみます．
>>> sp.std(sample_mean_array, ddof = 1)
0.2513580850860402
母標準偏差は0.8なので，母標準偏差よりも小さくなっていることがわかります．
標本平均のヒストグラム（平均値の標本分布）を描いてみると，以下のようになります．
>>> sns.distplot(sample_mean_array, color = 'black')
<matplotlib.axes._subplots.AxesSubplot object at 0x11b7ea0d0>
サンプルサイズが大きくなると標本平均は母平均に近くなることの確認
サンプルサイズと標本平均と母平均の関係をみてみます．シミュレーションは以下の条件で行います．
- 対象：標本平均
- 変化させるもの：サンプルサイズ
- 調べること：サンプルサイズが大きくなった時に，標本平均と母平均の関係はどうなるか
まず，10 〜 100010まで100区切りで変化させたサンプルサイズを用意します．
>>> size_array = np.arange(start = 10, stop = 100100, step = 100)
>>> size_array
array([ 10, 110, 210, ..., 99810, 99910, 100010])
標本平均を格納する配列を準備します．
>>> sample_mean_array_size = np.zeros(len(size_array))
シミュレーションを実行します（実行は，標本平均を求める試行を，サンプルサイズを変えながら何度も実行する）．
>>> np.random.seed(1)
>>> for i in range(0, len(size_array)):
... sample = population.rvs(size = size_array[i])
... sample_mean_array_size[i] = sp.mean(sample)
横軸にサンプルサイズをとって，縦軸に標本平均をとった折れ線グラフを描いてみます．
>>> plt.plot(size_array, sample_mean_array_size, color = 'black')
[<matplotlib.lines.Line2D object at 0x11cdd4b10>]
>>> plt.xlabel("sample size")
Text(0.5, 0, 'sample size')
>>> plt.ylabel("sample mean")
Text(0, 0.5, 'sample mean')
このグラフから，サンプルサイズが大きくなるにつれて，標本平均が母平均（4）に近づくことがわかります．
標本平均を何度も計算する関数
シミュレーションを行うために，標本平均を何度も計算する関数（calc_sample_mean()関数）を作ってみます．
>>> def calc_sample_mean(size, n_trial):
... sample_mean_array = np.zeros(n_trial)
... for i in range(0, n_trial):
... sample = population.rvs(size = size)
... sample_mean_array[i] = sp.mean(sample)
... return(sample_mean_array)
この関数の特徴は，引数としてサンプルサイズ（size）と試行回数（n_trial）を指定できるようにしてあることです．このようにしておけば，サンプルサイズや試行回数を変えて，標本平均を何度も計算することが可能になります．なお，標本平均は試行回数の数（n_trial個）だけ得られることになります．
上記のコードの説明：
2行目：試行回数分の標本平均を格納する配列を準備しています．
3行目：繰り返し構文．試行回数だけ繰り返すようにします．
4行目：平均 4，標準偏差 0.8 の正規分布に従う母集団から，size個の標本を抽出します．
5行目：sample_mean_array に標本平均を格納します．
6行目：n_trial回計算された標本平均を返します．
サンプルサイズを変化させた時の標本平均の分布
動作確認（データを10個選んで標本平均を得る試行を10000回繰り返す）
>>> np.random.seed(1)
>>> sp.mean(calc_sample_mean(size = 10, n_trial = 10000))
4.004202422791747
この関数を用いて，サンプルサイズの影響をより詳細に確かめてみます．
確かめる方法は，サンプルサイズを10, 20, 30と変化させた時の標本平均の分布をバイオリンプロットで確認することとします．試行回数は10000回とします．
乱数の種を指定
サンプルサイズ 10 の場合
>>>
>>> size_10 = calc_sample_mean(size = 10, n_trial = 10000)
>>> size_10_df = pd.DataFrame({
... "sample_mean":size_10,
... "size" :np.tile("size 10", 10000)
... })
サンプルサイズ 20の場合
>>> size_20 = calc_sample_mean(size = 20, n_trial = 10000)
>>> size_20_df = pd.DataFrame({
... "sample_mean":size_20,
... "size" :np.tile("size 20", 10000)
... })
サンプルサイズ 30の場合
>>> size_30 = calc_sample_mean(size = 30, n_trial = 10000)
>>> size_30_df = pd.DataFrame({
... "sample_mean":size_30,
... "size" :np.tile("size 30", 10000)
... })
各結果を結合する．
>>> sim_result = pd.concat(
... [size_10_df, size_20_df, size_30_df])
結果の表示（最初だけ）
>>> print(sim_result.head())
sample_mean size
0 3.922287 size 10
1 3.864329 size 10
2 4.069530 size 10
3 3.857140 size 10
4 4.184654 size 10
上記のコードでは，サンプルサイズを10, 20, 30と変化させた時の標本平均をpandasフレームにまとめています．
また，各場合（10, 20, 30）の試行回数が10000回（10000行）なので，結合した sim_result は30000行になります．
このデータをバイオリンプロットで描いてみます．
>>> sns.violinplot(x = "size", y = "sample_mean", data = sim_result, color = 'gray')
この図から，サンプルサイズが大きくなると，標本平均のばらつきが小さくなり，母平均（4）の近くに集中する様子がわかります．
標本平均の標準偏差が母標準偏差よりも小さいこと
サンプルサイズが大きくなると標本平均のばらつきが小さくなることを，標本平均の標準偏差をサンプルサイズ別にみてみることで確認してみます．シミュレーションは以下の要領で行います．
- 対象：標本平均の標準偏差
- 変化させるもの：サンプルサイズ
- 調べる内容：サンプルサイズが大きくなると，標本平均の標準偏差は小さくなっていくこと．すなわち，サンプルサイズを大きくすると，標本平均はより信頼できる値になっていく．
まずは，2- 100まで2区切りで変化させたサンプルサイズを用意します．
>>> size_array = np.arange(
... start =2, stop = 102, step =2)
>>> size_array
array([ 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26,
28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52,
54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78,
80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100])
続いて，標本平均の標準偏差を格納する入れ物を用意します．
>>> sample_mean_std_array = np.zeros(len(size_array))
シミュレーションを，試行回数は100回とします．
>>> np.random.seed(1)
>>> for i in range(0, len(size_array)):
... sample_mean = calc_sample_mean(size = size_array[i], n_trial = 100)
... sample_mean_std_array[i] = sp.std(sample_mean, ddof = 1)
横軸にサンプルサイズをとり，縦軸に標本平均の標準偏差をとった折れ線グラフを描いてみます．
>>> plt.plot(size_array, sample_mean_std_array, color = 'black')
[<matplotlib.lines.Line2D object at 0x119c87ed0>]
>>> plt.xlabel("sample size")
Text(0.5, 0, 'sample size')
>>> plt.ylabel("mean_std value")
Text(0, 0.5, 'mean_std value')
標準誤差
理論上の標本平均の標準偏差の大きさは，計算することができ，これを標準誤差（Standard Error, SE）と呼びます．
標準誤差（SE）$= \cfrac{\sigma}{\sqrt{N}}$
標準誤差とシミュレーションの結果を比較してみます．
まずは，標準誤差を求めます．
>>> standard_error = 0.8 / np.sqrt(size_array)
>>> standard_error
array([0.56568542, 0.4 , 0.32659863, 0.28284271, 0.25298221,
0.23094011, 0.21380899, 0.2 , 0.18856181, 0.17888544,
0.17056057, 0.16329932, 0.15689291, 0.15118579, 0.14605935,
0.14142136, 0.13719887, 0.13333333, 0.12977714, 0.12649111,
0.12344268, 0.12060454, 0.11795356, 0.11547005, 0.11313708,
0.11094004, 0.10886621, 0.1069045 , 0.10504515, 0.10327956,
0.1016001 , 0.1 , 0.09847319, 0.09701425, 0.09561829,
0.0942809 , 0.09299811, 0.09176629, 0.09058216, 0.08944272,
0.08834522, 0.08728716, 0.08626622, 0.08528029, 0.0843274 ,
0.08340577, 0.0825137 , 0.08164966, 0.0808122 , 0.08 ])
先の折れ線グラフに標準誤差を重ね合わせて描画してみます．
>>> plt.plot(size_array, sample_mean_std_array, color = 'black')
>>> plt.plot(size_array, standard_error, color = 'black', linestyle = 'dotted')
>>> plt.xlabel("sample size")
>>> plt.ylabel("mean_std value")
Text(0, 0.5, 'mean_std value')
この例では，linestyle = 'dotted' と点線で描くように指定しています．シミュレーション結果と標準誤差の値は，ほぼ一致しています．
標本分散の平均値と母分散のずれ
標本分散を対象としたシミュレーション（10000回）を実行して，標本分散の平均値を求めてみます．
まずは，標本分散を格納する変数（入れ物）を準備します．
>>> sample_var_array = np.zeros(10000)
シミュレーションを実行します．シミュレーションは，データを10個選んで，標本分散を求める思考としています．
>>> np.random.seed(1)
>>> for i in range(0, 10000):
... sample = population.rvs(size = 10)
... sample_var_array[i] = sp.var(sample, ddof = 0)
標本分散の平均値は以下のようになります．
>>> sp.mean(sample_var_array)
0.5746886877332101
母分散は0.8の2乗で0.64となるはずですが，計算結果は0.575と分散を過小評価していることがわかります．
不偏分散を使うとバイアスが解消される
標本分散には，母集団の分散（母分散）と比べて分散を過小評価してしまうというバイアスがあります．不偏分散は，これを修正するためのものです．
以下に，不偏分散を使ってそのことを確認してみます．
不偏分散を格納する変数（入れ物）を準備する．
>>> unbias_var_array = np.zeros(10000)
データを10個選んで不偏分散を求める試行を10000回繰り返す．
>>> np.random.seed(1)
>>> for i in range(0, 10000):
... sample = population.rvs(size = 10)
... unbias_var_array[i] = sp.var(sample, ddof = 1)
不偏分散の平均値を求める．
>>> sp.mean(unbias_var_array)
0.6385429863702334
先の標本分散を用いた結果よりも，不偏分散を用いた結果の方が，0.64に近い値となっていることがわかります．
サンプルサイズが大きくなると不偏分散は母分散に近づく
サンプルサイズと母分散の関係を調べるためのシミュレーションを行ってみます．
シミュレーションは以下の要領で行ってみます．
- 対象：不偏分散
- 変化させるもの：サンプルサイズ
- 調べる内容：サンプルサイズが大きくなると不偏分散は母分散に近づくこと
まずは，10 - 100010 まで100区切りで変化させたサンプルサイズを用意します．
>>> size_array = np.arange(
... start = 10, stop = 100100, step = 100)
>>> size_array
array([ 10, 110, 210, ..., 99810, 99910, 100010])
続いて不偏分散を格納する変数（入れ物）を用意します．
>>> unbias_var_array_size = np.zeros(len(size_array))
不偏分散を求める試行をサンプルサイズを変えて実行します．
>>> np.random.seed(1)
>>> for i in range(0, len(size_array)):
... sample = population.rvs(size = size_array[i])
... unbias_var_array_size[i] = sp.var(sample, ddof = 1)
横軸にサンプルサイズをとって，縦軸に不偏分散をとった折れ線グラフを描いてみます．
>>> plt.plot(size_array, unbias_var_array_size, color = 'black')
[<matplotlib.lines.Line2D object at 0x119c95fd0>]
>>> plt.xlabel("sample size")
Text(0.5, 0, 'sample size')
>>> plt.ylabel("unbias var")
Text(0, 0.5, 'unbias var')
この図から，サンプルサイズが大きくなるにつれて，不偏分散が母分散（0.64）に近ずく様子がわかります．
不偏性：推定量の期待値が真の母数（母集団のパラメタ）となる特性のこと．
不偏性があるということは，平均すると，過大にも過少にもなっていない，すなわち偏りがない推定量であるということになります．
一致性：サンプルサイズが大きくなると，推定量が真の母数に近づいていく特性
一致性があるということは，言い換えると，サンプルサイズが無限だった場合には，推定量と母数が一致するということになります．
分析の目的は，母集団分布を推定することです．母集団分布が分かれば，未知のデータに対して予測や推測を行うことができるからです．
母集団分布に正規分布を仮定すると，正規分布の母数（パラメータ）を推定することで，母集団分布を推定することができます．正規分布の母数は平均と分散の2つになります．正規分布の母数は一般的にはわからないので，代わりに標本平均と不偏分散を用いることになります．
上記のシミュレーションで，
- 標本平均と不偏分散が母数の推定量として好ましい性質を持つこと
がわかりました．また，
- 標本平均の平均値は母平均とほぼ同じであること
- 不偏分散の平均値は母分散とほぼ同じとみなせること
もわかりました．標本平均と不偏分散は不偏性を持っていると言うことができます．さらに，
- サンプルサイズを大きくすると，標本平均は母平均に近づくこと
- サンプルサイズを大きくすると，不偏分散は母分散に近づくこと
もわかりました．標本平均と不偏分散は一致性を持っていると言うことができます．
大数の法則：標本の大きさが大きくなるにつれて，標本平均が母平均に近づく近づき方を表現した法則，Upton, Cook (2010)
大数の法則には，大数の弱法則と大数の強法則がある．
中心極限定理：母集団分布が何であっても，サンプルサイズが大きい場合には，確率変数の和は正規分布に近いものになること．
中心極限定理をシミュレーションで確認してみます．シミュレーションは10000枚のコインを同時に投げて，表が出た枚数を数える試行を指定した回数だけ繰り返すことにします．
サンプルサイズと試行回数を指定します．
>>> n_size = 10000
>>> n_trial = 50000
表ならば1，裏ならば0とします．
>>> coin = np.array([0, 1])
表が出た回数
>>> count_coin = np.zeros(n_trial)
コインをn_size回投げる試行をn_trial回行います．
>>> np.random.seed(1)
>>> for i in range(0, n_trial):
... count_coin[i] = sp.sum(
... np.random.choice(coin, size = n_size, replace = True))
ヒストグラムを描いてみます．
>>> sns.distplot(count_coin, color = 'black')
<matplotlib.axes._subplots.AxesSubplot object at 0x115f41a10>
上記のコードをJupyter Notebookで実行した結果は
GitHubで見ることができます．
Pythonで学ぶあたらしい統計学の教科書に関連した過去の投稿
I 統計学の基礎 1. 統計学，2. 標本が得られるプロセス，3. 標本が得られるプロセスの抽象化
I 統計学の基礎 4. 記述統計の基礎，6. 確率質量関数と確率密度関数
I 統計学の基礎 7. 統計量の基礎
I 統計学の基礎 8. 確率論の基本
I 統計学の基礎 9. 確率変数と確率分布
III データ分析 記述統計：1変量データ編
III データ分析 記述統計：多変量データ編
III データ分析 母集団からの標本抽出シミュレーション
III 標本の統計量の性質
のまとめやコードを実行してみた際のメモです．
