中心極限定理のシミュレーション

 

中心極限定理によれば、ある程度の標本サイズがあれば抽出元の母集団が正規分布でないとしても標本平均は正規分布に従う。 本当だろうか。証明は難しいのでPythonでシミュレーションを行って確認する。

まずは与えられたデータのヒストグラム、要約、歪度、尖度、正規Q-Qプロットを表示する関数を定義する。

In [146]:
%matplotlib inline
from numpy.random import normal, poisson
import scipy.stats as stats

def show(data):
    data = Series(data)
    plt.hist(data)
    plt.show()
    print data.describe()
    print "skewness:", data.skew()
    print "kurtosis:", data.kurtosis()
    stats.probplot(data, dist="norm", plot=plt)
    plt.show()
 

平均16、標準偏差16の正規分布に従うランダムなデータを作り、この関数を適用してみる。

In [147]:
data = normal(loc=16, scale=16, size=10000)
show(data)
 
 
count    10000.000000
mean        15.959868
std         15.982171
min        -44.150475
25%          5.131434
50%         15.923589
75%         26.594607
max         76.952743
dtype: float64
skewness: 0.003154636326
kurtosis: 0.0415231202693
 
 

meanとstdから、この集団の平均と標準偏差が大体16であることがわかる。 歪度skewnessはゼロに近いほど左右対称であることを示す。尖度kurtosisもゼロに近いほど正規分布に近くなる。 最後の正規Q-Qプロットは青いデータが赤い線に乗っているほど正規分布に従っているということになる。

次に、データからサンプリングを行ってその平均の標本分布を表示する関数を定義する。

In [148]:
def sample_and_show_means(data, num_samples=100, num_iterations=1000):
    print "Sample Distribution"
    population = Series(data)
    sample_means = []
    for i in range(0, num_iterations):
        sample_means.append(population.sample(n=num_samples).mean())
    sample_means = Series(sample_means)
    show(sample_means)
 

この関数はデフォルトでは標本サイズ100のサンプリングを1000回行う。これを使って先ほどの正規分布の母集団からサンプリングを行う。

In [149]:
sample_and_show_means(data)
 
Sample Distribution
 
 
count    1000.000000
mean       15.928242
std         1.629540
min        11.022722
25%        14.780568
50%        15.910704
75%        17.111628
max        20.185614
dtype: float64
skewness: -0.0608301712331
kurtosis: -0.183327796726
 
 

このように標本平均の分布も正規分布に従っているようだ。標本平均の標準偏差は標準誤差と呼ばれる統計量で、標準偏差/(標本サイズの平方根)であるとされる。理論上は16/sqrt(100)=1.6ということになり、それに近い値が出ている。

標本サイズが小さくなるほど標準誤差は大きくなる。今度は標本サイズを16にしてみよう。

In [150]:
sample_and_show_means(data, num_samples=16)
 
Sample Distribution
 
 
count    1000.000000
mean       15.787418
std         3.971251
min         3.899584
25%        12.963414
50%        15.981024
75%        18.451290
max        28.748238
dtype: float64
skewness: -0.0268117549844
kurtosis: -0.118104482425
 
 

理論上の標準誤差は16/sqrt(16)=4で、やはりそれに近い値が出ている。

今度は正規分布でない母集団からサンプリングした場合も標本平均の分布が正規分布に従うのかどうかを確かめる。まずはポワソン分布を見る。

In [151]:
data = poisson(lam=2, size=10000)
show(data)
sample_and_show_means(data)
 
 
count    10000.000000
mean         2.010800
std          1.413578
min          0.000000
25%          1.000000
50%          2.000000
75%          3.000000
max         10.000000
dtype: float64
skewness: 0.709456227392
kurtosis: 0.594243907564
 
 
Sample Distribution
 
 
count    1000.000000
mean        2.023040
std         0.139944
min         1.590000
25%         1.930000
50%         2.020000
75%         2.120000
max         2.510000
dtype: float64
skewness: 0.0248246692832
kurtosis: 0.00103792770786
 
 

ポワソン分布自体は正規分布から遠いのに対して、標本平均の分布はかなり正確に正規分布に従っていることがわかる。

正規分布やポワソン分布のような由緒正しい分布ではなく、名前のつかないようなでたらめな分布でも同じなのだろうか。 今度は2つの分布を混ぜ合わせたものを母集団としてみる。

In [152]:
data = np.concatenate([poisson(lam=2, size=10000), normal(loc=16, scale=2, size=10000)])
show(data)
sample_and_show_means(data)
 
 
count    20000.000000
mean         9.002401
std          7.203268
min          0.000000
25%          2.000000
50%          9.000000
75%         15.967441
max         23.631625
dtype: float64
skewness: 0.0588665817353
kurtosis: -1.78155778418
 
 
Sample Distribution
 
 
count    1000.000000
mean        8.973318
std         0.732240
min         6.892384
25%         8.472600
50%         8.975025
75%         9.484283
max        11.293600
dtype: float64
skewness: -0.0226618272725
kurtosis: -0.220939870543
 
 

母集団のほうはなんだかよくわからない分布であるが、標本平均のほうはとにかく正規分布に従っていることがわかる。