Pythonで大数の法則からリンデブルグの中心極限定理まで

今回は中心極限定理を解説していきます。

中心極限定理は、大数の法則、チェビシェフの不等式を先に学ぶと理解しやすくなります。

大数の法則

これは、「標本数が多いほど、標本平均は母平均に近い値になる」ということを表す法則です。定理は以下の通りです。

定理：（大数の法則）互いに独立で同一の平均値$\mu$をもつ確率変数$X_1, X_2, ... X_n, ...$に対し、それぞれの分散が$ \sigma_j^2 \leq \sigma < \infty (j = 1, 2, ...)$を満たす（分散は異なっていてもよい）とき、以下の式が成り立つ

$$ \lim_{n → \infty} P(|\frac{X_1 + X_2 + ... + X_n}{n} - \mu | > \epsilon) = 0 $$

この$\epsilon$は誤差であり、$\epsilon = 1$などと任意の値を設定し、サンプルサイズ$n$を大きくしていった場合に、$[標本平均] - [母平均]$の絶対値が$\epsilon$よりも大きくなる確率が$0$に近づく、ということを表しています。

大数の法則が見られる面白い例として、モンテカルロ法や、大数の法則の暗号読解への応用があります。

モンテカルロ法は面積の計算法の一つで、$2 \cdot 2$の正方形に内接する円の面積を求める場合に、正方形内をランダムに選んで点を打つという試行を十分に繰り返すと、円の内側の点の数が全ての点の数の$\frac{\pi}{4}$に近づくため、円の面積が$\pi$になる、というものです。

Pythonで確かめてみます。

import numpy as np

# -1から1までの一様乱数x、yを作成
n = 10000
np.random.seed(seed=12345)
x = np.random.uniform(-1, 1, n)
y = np.random.uniform(-1, 1, n)

# 中心(0,0)、半径1の円内に入った点の数をカウント
p_in = 0
for k in range(n):
  if x[k]**2 + y[k]**2 < 1:
    p_in = p_in + 1

# 円内に入るのはπ/4の確率なので、4倍するとπになる
print((p_in / 10000) * 4)
# out: 3.1464

大数の法則の暗号読解への応用は、まずアルファベットを全て別のアルファベットに置き換えるという暗号を解読する場合、十分にサンプルがあれば、アルファベットeの出現回数がダントツで多く、次にtが多い、などとわかるため、すぐに解読できてしまうということです。これを頻度分析といいます。つまり統計学的な特徴が現れる暗号は解読されてしまうということです。

そしてこの大数の法則を定量的にあらわしたのが次のチェビシェフの不等式です。

チェビシェフの不等式

定理：（チェビシェフの不等式）平均$\mu$、分散$\sigma^2$をもつ確率変数$X$は、次の不等式を満たす。

$$ P(|X - \mu | > \epsilon) \leq \frac{\sigma^2}{\epsilon^2} $$

これは$\epsilon = k\sigma$とおいて、

$$ P(|X - \mu | > k\sigma) \leq \frac{1}{k^2} $$

と表記されることもあります。

$[0, 12]$を台にもつ一様乱数に対し、このチェビシェフの不等式を当てはめてみます。一様乱数$X$を$10000$個発生させて、標本平均を$Z_{30}$として、真の平均とのズレが$\epsilon = 1$より大きくなる確率を求めてみます。平均は$\mu = (0 + 12) / 2 = 6$、分散は$\sigma^2 = (12 - 0)^2 / 12 = 12$なので、

$$ P(|Z_{30} - 6| > 1) \leq \frac{12}{30 \cdot 1^2} = 0.4 $$

となり、チェビシェフの不等式では、真の平均とのズレが$1$より大きくなる確率は$0.4$以下になることがわかりました。Pythonでデータを作り、実測値としての確率を求めてみます。

np.random.seed(seed=12345)

# 10000行50列の[0,12]を台にもつ一様乱数を作成
x = np.random.uniform(0, 12, 30*10000)
xm = x.reshape(10000,30)

# 一行ごとの平均をz30に格納
def l_mean(x_i):
  return np.mean(x_i)
  
z30 = np.apply_along_axis(l_mean, 1, xm)

#標本平均と真の平均とのズレが1より大きくなる確率
print(sum(np.abs(z30 - 6)>1)/10000)
# out: 0.1118

このように、理論値$0.4$と実測値$0.1118$は差分が大きく、チェビシェフの不等式は真ではあるものの、精度があまりよくありません。この問題を解決したのが中心極限定理です。

中心極限定理

定理：（中心極限定理）平均$\mu$、分散$\sigma^2$の独立同分布の確率変数$X_1, X_2, ... X_n$に対して、$S_n = \sum_{j=1}^{n} X_j$とすると、以下が成り立つ

$$ P(a \leq \frac{S_n - n\mu}{\sqrt{n}\sigma} \leq b) = \int_{a}^{b} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}dx $$

まず連続確率変数は確率密度関数$f(x)$を用いて以下のように表すことができます。

$$ P(a \leq X \leq b) = \int_{a}^{b}f(x)dx $$

中心極限定理の式の右辺$\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$の部分を見ると、これは正規分布$\int_{a}^{b} \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x - \mu)^2}{2\sigma^2}}dx$の平均$\mu$が$0$、分散$\sigma^2$が$1$の正規分布、つまり標準正規分布$N(0, 1)$であることがわかります。

よって、$Y_n = \frac{S_n - n\mu}{\sqrt{n}\sigma}$の積率母関数が$N(0, 1)$の積率母関数に収束すれば、中心極限定理を証明できることになります。証明の過程は省略します。

中心極限定理は簡単に言えば、サンプルサイズ$n$が十分に大きい時、母集団の分布がどのようなものであっても、和$S_n = X_1 + X_2 + ... + X_n$の分布が$N(n\mu, n\sigma^2)$に近似し、標本平均$S_n/n$の分布が$N(\mu, \frac{\sigma^2}{n})$に近似するということです。

Pythonで指数分布の場合で考えてみます。まず指数分布に従うサンプルを300000個抽出します。

import matplotlib.pyplot as plt
# ある事象がある期間内に起きる回数の平均が10回とする
target = 10
beta = 1.0/target
np.random.seed(seed=12345)
x = np.random.exponential(beta, 10000*30)

plt.hist(x, density=True, bins=100,lw=0,alpha=.8)
plt.ylim(0,10)

この指数分布に従う確率変数$X_i (i=1,2,...10000)$を作成します。pandasのDataFrameで表示すると以下のようになります。

import pandas as pd
# 10000行30列に変形
xm = x.reshape(10000, 30)
pd.DataFrame(xm)

$Y_n = \frac{S_n - n\mu}{\sqrt{n}\sigma}$を計算し、標準正規分布$N(0,1)$と比較します。

from scipy.stats import norm

def f_Y(x_i):
  s = np.sum(x_i)
  n = len(x_i)
  mu = 1/target
  sigma = np.sqrt(1/target**2)
  return (s - n * mu) / (np.sqrt(n) * sigma)
  
y30 = np.apply_along_axis(f_Y, 1, xm)

# pdfで確率密度関数を生成
X = np.arange(start=-3, stop=3, step=0.001)
norm_pdf = norm.pdf(x=X, loc=0, scale=1)

# 可視化
plt.plot(X, norm_pdf)
plt.hist(y30, density=True, bins=100,lw=0,alpha=.8)

print(np.mean(y30))
print(np.var(y30))
# out: 0.00784183009181524
# out: 0.9831708609225359

$Y_n = \frac{S_n - n\mu}{\sqrt{n}\sigma}$の実測値の平均が$0$、分散が$1$に近似し、標準正規分布$N(0, 1)$に近似しているのがわかります。

つぎに平均の分布と、確率密度関数を重ねて見てみます。

def l_mean(X_i):
  return np.mean(X_i)
  
z30 = np.apply_along_axis(l_mean, 1, xm)

# pdfで確率密度関数を生成
X = np.arange(start=0, stop=0.2, step=0.001)
norm_pdf = norm.pdf(x=X, loc=np.mean(z30), scale=np.std(z30))

# 可視化
plt.plot(X, norm_pdf)
plt.hist(z30, density=True, bins=100,lw=0,alpha=.8)

print(np.mean(z30))
print(np.var(z30))
# out: 0.100143171574447
# out: 0.00032772362030751204

指数分布の期待値は$E(X) = \frac{1}{\lambda}$、分散は$V(X) = \frac{1}{\lambda^2}$です。今回は$\lambda = 10$に設定しています。$\frac{\sigma^2}{n}$は$\frac{0.01}{30} = 0.000333$となり、標本平均$S_n/n$の分布が$N(\mu, \frac{\sigma^2}{n})$に近似していることが確認できます。

中心極限定理の驚くべき点は、期待値と分散が存在しさえすれば、どのような分布でも、標本平均の分布がサンプルサイズを大きくするほど正規分布に近づくことです。

中心極限定理は強力ですが、同分布を仮定しています。中心極限定理を同分布を仮定しない場合に成り立つ形にした定理があります。それがリンデブルグの中心極限定理です。