[Statistics] III 母集団からの標本抽出シミュレーション

Pythonで学ぶあたらしい統計学の教科書のまとめやコードを実行してみた際のメモです．

まずは，数値計算に必要なライブラリをインポートします．

>>> 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

これもREPLでは必要ありませんが，Jupyter Notebookでグラフを表示させる場合には以下のようにします．

>>> %matplotlib inline

例として，5つしかデータがない場合を考えます．以下は，5つの母集団の設定です．

>>> fish_5 = np.array([2, 3, 4, 5, 6])

>>> fish_5

array([2, 3, 4, 5, 6])

母集団の中からランダムに抽出を行うには，np.random.choice 関数を用います．size = 1 とすることでサンプルサイズが1となります．replace = False と指定することで，同じ標本が2回選ばれることがない設定としています．

>>> np.random.choice(fish_5, size = 1, replace = False)

array([2])

続いて，3つの値をサンプルリングしてみます．

>>> np.random.choice(fish_5, size = 3, replace = False)

array([3, 6, 4])

もう一度サンプリングしてみます．

>>> np.random.choice(fish_5, size = 3, replace = False)

array([3, 2, 4])

確率的に得られる標本が変わるため，全く同じコード np.random.choice(fish_5, size = 3, replace = False を実行しても，毎回異なる標本が得られています．np.random.choice関数を使うと，毎回ランダムに標本が選ばれますが，乱数の種類を指定すると，毎回，同じデータがランダムに選ばれることになります．乱数の種類を指定するには，np.random.seed関数を用います．

以下の例では，乱数の種類に1（np,random.seed(1)）を指定して実行しています．

>>> np.random.seed(1)

>>> np.random.choice(fish_5, size = 3, replace = False)

array([4, 3, 6])

>>> np.random.seed(1)

>>> np.random.choice(fish_5, size = 3, replace = False)

array([4, 3, 6])

先の例と違って，全く同じ標本が得られていますが，これは乱数の種を指定した効果です．

得られた標本から平均値を計算して，標本平均を求めてみます．

>>> np.random.seed(1)

>>> sp.mean(np.random.choice(fish_5, size = 3, replace = False))

4.333333333333333

乱数：ランダムに得られる値のこと（分野によっては確率変数とほぼ同様の意味で用いられる）

復元抽出：抽出した標本をまた母集団に戻してから再度抽出するサンプリングの方法のこと．

非復元抽出：抽出した標本を元に戻すことはしないサンプリング方法のこと．

上記のコードでは，replace = False の部分が復元抽出なのか，非復元抽出なのかを指定している部分になります．

ここまでの例は，母集団の数が5つの例でしたが，ここからは母集団の数を10000にしてみます（10000個のデータの入っている3-4-1-fish_length_100000.csv* ファイルを読み込みます）．

>>> fish_10000 = pd.read_csv("3-4-1-fish_length_100000.csv")["length"]

>>> fish_10000.head()

0    5.297442

1    3.505566

2    3.572546

3    3.135979

4    4.689275

Name: length, dtype: float64

データ（母集団）の個数を確認します．

>>> len(fish_100000)

100000

母集団のデータの個数は増えましたが，標本を抽出する方法は，先のものと同じです．

>>> sampling_result = np.random.choice(fish_100000, size = 10, replace = False)

>>> sampling_result

array([3.35462508, 4.32985663, 4.10081464, 4.0318259 , 3.80132081,

       4.37080993, 3.78355408, 6.59080065, 4.90904882, 3.54286598])

標本平均を計算してみます．

>>> sp.mean(sampling_result)

4.281552250352997

母集団の数が増えたので，その平均値（母平均）を求めてみます．

>>> sp.mean(fish_100000)

4.000000000000001

母標準偏差を求めてみます．

>>> sp.std(fish_100000, ddof = 0)

0.8

母分散を求めてみます．

>>> sp.var(fish_100000, ddof = 0)

0.64

母集団のヒストグラムを描いてみると以下のようになります．

>>> sns.distplot(fish_100000, kde = False, color = 'black')

ヒストグラムを見ると，平均値（4）を中心として，左右対称な度数分布になっているのがわかります．ここで，母集団の確率分布は，平均4.00，分散0.64の正規分布として表現できると仮定して，以下に母集団分布と正規分布の確率密度関数の比較を行ってみます．

平均4.00，分散0.64の正規分布の確率密度を，1〜7の範囲で図示してみます．

まずは，1〜7までを0.1区切りで分けた等差数列を用意します．

>>> x = np.arange(start = 1, stop = 7.1, step = 0.1)

>>> x

array([1. , 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2,

       2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3. , 3.1, 3.2, 3.3, 3.4, 3.5,

       3.6, 3.7, 3.8, 3.9, 4. , 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8,

       4.9, 5. , 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6. , 6.1,

       6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7. ])

確率密度は stats.norm.pdf 関数を用いて計算できます．引数 loc が平均値で，scaleが標準偏差です．

>>> stats.norm.pdf(x = x, loc = 4, scale = 0.8)

array([4.40744603e-04, 6.98826903e-04, 1.09085337e-03, 1.67639859e-03,

       2.53631007e-03, 3.77782254e-03, 5.53981051e-03, 7.99765039e-03,

       1.13669531e-02, 1.59052270e-02, 2.19103756e-02, 2.97148760e-02,

       3.96745648e-02, 5.21512316e-02, 6.74887081e-02, 8.59828448e-02,

       1.07846649e-01, 1.33172835e-01, 1.61896995e-01, 1.93765332e-01,

       2.28311357e-01, 2.64845807e-01, 3.02463406e-01, 3.40068748e-01,

       3.76421790e-01, 4.10201211e-01, 4.40081658e-01, 4.64818867e-01,

       4.83335146e-01, 4.94797109e-01, 4.98677851e-01, 4.94797109e-01,

       4.83335146e-01, 4.64818867e-01, 4.40081658e-01, 4.10201211e-01,

       3.76421790e-01, 3.40068748e-01, 3.02463406e-01, 2.64845807e-01,

       2.28311357e-01, 1.93765332e-01, 1.61896995e-01, 1.33172835e-01,

       1.07846649e-01, 8.59828448e-02, 6.74887081e-02, 5.21512316e-02,

       3.96745648e-02, 2.97148760e-02, 2.19103756e-02, 1.59052270e-02,

       1.13669531e-02, 7.99765039e-03, 5.53981051e-03, 3.77782254e-03,

       2.53631007e-03, 1.67639859e-03, 1.09085337e-03, 6.98826903e-04,

       4.40744603e-04])

>>> plt.plot(x, stats.norm.pdf(x = x, loc = 4, scale = 0.8), color = 'black')

[<matplotlib.lines.Line2D object at 0x10edef048>]

上記で表示されている配列における e-04 は，10のマイナス4乗を意味しています．

続いて，正規分布の確率密度と母集団のヒストグラムのグラフを重ねてみます．

>>> sns.distplot(fish_100000, kde = False, norm_hist = True, color = 'black')

<matplotlib.axes._subplots.AxesSubplot object at 0x12bd9d588>

>>> plt.plot(x, stats.norm.pdf(x = x, loc = 4, scale = 0.8), color = 'black')

[<matplotlib.lines.Line2D object at 0x11bfc2a58>]

上記では，sns.distplot において，norm_hist = True と指定することで，面積が1であるヒストグラムになるように指定しています．

上図から，母集団分布は平均4.00, 分散0.64の正規分布とみなすことができそうだと推定できます．このことは，母集団からの標本抽出は，正規分布に従う乱数を発生させることとほぼ同じだとみなすことができるということを意味します．

今回の投稿の例では，最初は母集団として fish_100000 を用いてそのデータからnp.random.choice 関数を使って標本抽出していましたが，それらを使わずに，最初から正規分布に従う乱数を発生させる関数（stats.norm.rvs 関数）を使用してみます．
stats.norm.rvs 関数では，平均 loc, 標準偏差 scale, 取得するサンプルサイズ size の3つを引数として取得します．

>>> sampling_norm = stats.norm.rvs(loc = 4, scale = 0.8, size = 10)

>>> sampling_norm

array([4.53940856, 4.79027469, 4.35494541, 4.61073595, 3.27863717,

       4.35708691, 3.22347338, 2.72309008, 3.99801584, 4.26737518])

標本平均を求めることもできます．

>>> sp.mean(sampling_norm)

4.014304314985844

上記の例で，stats.norm.rvs 関数を用いて正規乱数を生成することは，母集団からの標本抽出を行っていることと同じことになります．

母集団のヒストグラムと正規分布の確率密度が等しいとみなしたことについて
ヒストグラムは階級ごとに度数を図示しているので，当然ながらグラフはギザギザの形になります．一方で，正規分布の密度関数は滑らかです．よって，厳密には正規分布の確率密度とヒストグラムは完全に一致していません．
母集団分布に正規分布を仮定するということは，母集団が無限母集団であり，級数を無限に増やすことで，無限に細かいヒストグラムを描いたとすると正規分布の確率密度関数と一致するということになります．

これらのことをJupyter Notebookで実行した結果はGitHubで見ることができます．