Python
確率分布
統計学
中心極限定理


今回のブログの趣旨


  • 統計検定2級レベルで頻出の確率分布について改めてメモとしてまとめる!

  • 実際にシミュレーションしてみることで、中心極限定理を実感する!


中心極限定理とは

母集団の平均を $\mu$、分散を$\sigma^2$ とした時、母集団の分布がどんな分布であっても、そこから $n$ 個ランダムサンプリングした標本の平均は、$n$ が十分に大きい時、近似的に平均 $\mu$、分散 $\sigma^2/n$ の正規分布に従います。尚、母平均と標本平均の誤差を標準誤差(SE:Standard Error)といいます。


様々な確率分布

今回は以下の確率分布において、中心極限定理をシミュレーションしてみます。

関数
分布
連続/離散

uniform
一様分布
連続

binomial
二項分布
連続

poisson
ポアソン分布
連続

まず、必要なライブラリと変数を宣言します。

import scipy as sp

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
sns.set()
%matplotlib inline


一様分布


  • いかなる確率変数でも確率密度関数の値が、等しい分布

  • 以下に0から100までの一様分布(離散値)を図示してみます。


x = np.linspace(0, 100, 10000)
mean = x.mean()
mean_line_x = np.array([mean, mean])
mean_line_y = np.array([0, 1000])

# 一様分布のヒストグラムを可視化
plt.hist(x, alpha=0.5)

# 分布の平均と分散を表示
print('x_mean: ', x.mean(), '\nx_var: ', x.var())

# 分布の平均をプロット
plt.plot(mean_line_x, mean_line_y, color='r', linewidth=2, linestyle = "--", label = "mean")
plt.legend();

download.png

因みに赤い線は分布の平均(期待値)です。この時の平均値と分散は以下の通りでした。

$\mu=50.00$

$\sigma^2=833.50$

次にこの一様分布から、「10個のサンプルを取り平均を求める」という作業を1万回やった時のサンプル平均の分布を図示化してみます。

N = 10000

n = 10

np.random.seed(0)
x = np.array([])

for _ in range(N):
# 0から100までの一様分布(連続値)からサンプル数nを抽出
uniform = np.random.uniform(0, 100, n)
sample_mean = uniform.mean()
x = np.append(x, sample_mean)

mean_x = x.mean()
mean_line_x = np.array([mean, mean])
mean_line_y = np.array([0, 3000])

# 分布の平均と分散を表示
print('x_mean: ', x.mean(), '\nx_var: ', x.var(ddof=1))

# 標本平均の分布とその平均を表示
plt.hist(x, alpha=0.5)
plt.plot(mean_line_x, mean_line_y, color='r', linewidth=2, linestyle = "--", label = "mean")
plt.legend();

download-1.png

見事に釣鐘型になりました。因みにこの分布の平均と不偏分散は以下の通りで、「標本の不偏分散は $\sigma^2/n$ (母分散/標本のサンプル数)になる。」というのも何となく確認出来ました。

$\mu=49.95$

$\sigma^2=83.80$


二項分布

「コインを投げたときに表が出るか裏が出るか」のように、何かを行ったときに起こる結果が2つしかない試行のことを「ベルヌーイ試行」といいます。このベルヌーイ試行を数回行って、成功する回数が従う確率分布を「二項分布」といいます。今、試行回数を $n$、確率変数(成功する回数)を $X$、成功確率を $p$とすると

$$ P(X) = nCx×p^x(1-p)^{n-x}$$

となります。

from scipy.special import comb

n = 100 # 試行回数
p = 0.3 # 成功確率

x = np.linspace(0, 100, 1)
xList = pd.Series(
[comb(float(n), x)*p**x*(1-p)**(float(n)-x) for x in range(0, n+1)])

mean = n * p
var = n * p * (1 - p)
print('mean: ',mean, '\nvar: ', var)

mean_line_x = np.array([mean, mean])
mean_line_y = np.array([0, 0.09])

plt.plot(xList.index, xList)
plt.plot(mean_line_x, mean_line_y, color='r', linewidth=2, linestyle = "--", label = "mean")
plt.legend();

download-2.png

$\mu=30.0$

$\sigma^2=21.0$


  • 二項分布の試行回数を大きくしていくと正規分布に近づくことも何となく確認できます。

  • 次にこの二項分布から、「10個のサンプルを取り平均を求める」という作業を1万回やった時のサンプル平均の分布を図示化してみます。

N = 10000

n = 10

np.random.seed(0)
x = np.array([])

for _ in range(N):
# 0から100までの一様分布(連続値)からサンプル数nを抽出
binomial = np.random.binomial(100, 0.3, n)
sample_mean = binomial.mean()
x = np.append(x, sample_mean)

mean_x = x.mean()
mean_line_x = np.array([mean, mean])
mean_line_y = np.array([0, 3000])

# 分布の平均と分散を表示
print('x_mean: ', x.mean(), '\nx_var: ', x.var(ddof=1))

# 標本平均の分布とその平均を表示
plt.hist(x, alpha=0.5)
plt.plot(mean_line_x, mean_line_y, color='r', linewidth=2, linestyle = "--", label = "mean")
plt.legend();

download-3.png

やはり中心極限定理を何となく確認できました。


ポアソン分布


  • ポアソン分布とは、「単位時間あたりに平均 λ 回起こる現象が、単位時間に k 回起きる確率」の分布となります。

  • 二項分布との違いはn(試行回数)が非常に大きく、逆にp(成功確率)が非常に小さいことです。

  • 平均 $ \lambda = np $

  • 分散 $ \lambda = np $

  • 平均と分散が等しくなるという性質を持っています

  • 二項分布の分散 $np(1-p)$ で $n$ を非常に大きく $p$ を非常に小さくすると $np$ となることからも上記は分かると思います。

  • 1年あたり平均0.61人の兵士が馬に蹴られて死ぬ軍隊において、「1年に何人の兵士が馬に蹴られて死ぬかの確率の分布」を求める。それが、歴史上で初めてポアソン分布が使われた事例だと言われています。(メモ)

$\lambda$=10のポアソン分布を図示します

import math

n = 100
lam = 10

x = np.linspace(0, n, 1)
xList = pd.Series(
[((math.e ** -lam) * (lam ** x)) / math.factorial(x) for x in range(0, n+1)])

mean = lam
var = lam
print('mean: ',mean, '\nvar: ', var)

mean_line_x = np.array([mean, mean])
mean_line_y = np.array([0, 0.13])

plt.plot(xList.index, xList)
plt.plot(mean_line_x, mean_line_y, color='r', linewidth=2, linestyle = "--", label = "mean")
plt.legend();

download-4.png

次に、このポアソン分布から、「10個のサンプルを取り平均を求める」という作業を1万回やった時のサンプル平均の分布を図示化してみます。

N = 10000

n = 10

np.random.seed(0)
x = np.array([])

for _ in range(N):
# 0から100までの一様分布(連続値)からサンプル数nを抽出
poisson = np.random.poisson(10, n)
sample_mean = poisson.mean()
x = np.append(x, sample_mean)

mean_x = x.mean()
mean_line_x = np.array([mean, mean])
mean_line_y = np.array([0, 3000])

# 分布の平均と分散を表示
print('x_mean: ', x.mean(), '\nx_var: ', x.var(ddof=1))

# 標本平均の分布とその平均を表示
plt.hist(x, alpha=0.5)
plt.plot(mean_line_x, mean_line_y, color='r', linewidth=2, linestyle = "--", label = "mean")
plt.legend();

download-5.png

先ほどまでより、少し分かりづらいですが、こちらも正規分布に近づいている様子がみて取れます。


今回のブログで学んだこと


  • 中心極限定理を確認した

  • matplotlibを手足のように使いこなすのは難しい