2024.06.01
数学

中心極限定理

はじめに

　今回は中心極限定理について説明する。数学的な証明は行わず、その振る舞いについてPythonプログラムで確認したあと、この定理の有用性の例を示す。

中心極限定理とは

　任意の確率分布 $p(x)$ を考える。その平均 $\mu$ と分散 $\sigma^2$ は既知であるとする。この確率分布から独立に $N$ 個のデータをサンプリングする。

$\begin{align*} D=\{x_1,\cdots,x_N\} \end{align*}$

「独立に」とは、今のサンプリングが次のサンプリングに影響を与えないということである。このような $N$ 個のサンプリングをたくさん繰り返す。

$\begin{align*} D^{(1)}=\{x_1^{(1)},\cdots,x_N^{(1)}\},\cdots, D^{(K)}=\{x_1^{(K)},\cdots,x_N^{(K)}\} \end{align*}$

各データの組 $D^{(k)}$ から平均（標本平均）を計算する。

$\begin{align*} \bar{x}^{(k)}=\frac{1}{N}\sum_{n=1}^N x_n^{(k)},\;\;\;k=1,\cdots,K \end{align*}$

この平均の分布 $\{\bar{x}^{(1)},\cdots,\bar{x}^{(K)}\}$ は、 $N$ を十分に大きくとると、平均が $\mu$ 、分散が $\sigma^2/N$ の正規分布に近づくことが証明されている。これが中心極限定理である。ここで大切なことは、最初に考えた確率分布 $p(x)$ の種類は任意であることである。ベルヌーイ分布でも二項分布でも構わない。サンプリング数を十分大きくすればその平均は正規分布に近づくのである（下図参照）。

図1

Pythonによる実験

　ここでは任意の確率分布として0以上1未満で定義された一様分布 $p(x)$ を考える。

$\begin{align*} p(x)=1 \end{align*}$

このとき

$\begin{align*} \int^1_0dx\;p(x)=1 \end{align*}$

が成り立ち、平均と分散は

$\begin{align*} \mu&=\int^1_0dx\;p(x)\;x = \frac{1}{2} \\ \sigma^2&=\int^1_0dx\;p(x)\;(x-\mu)^2 = \frac{1}{12} \end{align*}$

となる。さて、最初に以下の手順を実装する。

一様分布から $N$ 個のデータを無作為に取り出す。
$N$ 個のデータの平均を計算する。
1と2を $K$ 回繰り返す。

def calcuate_means(N: int, K: int) -> list[float]:
    x_means = []
    for _ in range(K):
        xs = []
        for n in range(N):
            x = np.random.rand()  # 一様分布からの乱数生成 [0,1)
            xs.append(x)
        mean = np.mean(xs)  # N個のサンプルの平均
        x_means.append(mean)  # K回繰り返す。
    return x_means

def calcuate_means(N: int, K: int) -> list[float]:

x_means = []

for _ in range(K):

xs = []

for n in range(N):

x = np.random.rand() # 一様分布からの乱数生成 [0,1)

xs.append(x)

mean = np.mean(xs) # N個のサンプルの平均

x_means.append(mean) # K回繰り返す。

return x_means

残りは、描画部分とメイン関数部分である。

def draw_graph(x_means: list[float], N: int, K: int, path: str, mu: float, sigma: float) -> None:
    plt.hist(x_means, bins="auto", density=True)
    plt.title(f"N={N}, K={K}")
    plt.xlabel("x")
    plt.ylabel("Probability Density")
    plt.xlim(-0.05, 1.05)
    plt.ylim(0, 5)
    # draw normal distribution
    x = np.linspace(-0.05, 1.05, 100)
    y = 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-((x - mu) ** 2) / (2 * sigma**2))
    plt.plot(x, y, color="red", linestyle="dashed")
    plt.savefig(path)
    plt.clf()


def make_normal(mu: float, sigma: float):
    x = np.random.normal(mu, sigma)
    return x


if __name__ == "__main__":
    K = 10000
    MU = 0.5
    SIGMA = np.sqrt(1.0 / 12.0)
    for n in [1, 2, 4, 10]:
        x_means = calcuate_means(n, K)
        path = f"outputs/pdf_{n:02d}.png"
        draw_graph(x_means, n, K, path, MU, SIGMA / np.sqrt(n))

def draw_graph(x_means: list[float], N: int, K: int, path: str, mu: float, sigma: float) -> None:

plt.hist(x_means, bins="auto", density=True)

plt.title(f"N={N}, K={K}")

plt.xlabel("x")

plt.ylabel("Probability Density")

plt.xlim(-0.05, 1.05)

plt.ylim(0, 5)

# draw normal distribution

x = np.linspace(-0.05, 1.05, 100)

y = 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-((x - mu) ** 2) / (2 * sigma**2))

plt.plot(x, y, color="red", linestyle="dashed")

plt.savefig(path)

plt.clf()

def make_normal(mu: float, sigma: float):

x = np.random.normal(mu, sigma)

return x

if __name__ == "__main__":

K = 10000

MU = 0.5

SIGMA = np.sqrt(1.0 / 12.0)

for n in [1, 2, 4, 10]:

x_means = calcuate_means(n, K)

path = f"outputs/pdf_{n:02d}.png"

draw_graph(x_means, n, K, path, MU, SIGMA / np.sqrt(n))

特に難しい部分はない。実験結果は以下の図の通り。図内の $N$ と $K$ は上の説明で使用したパラメータに対応する。 $K$ を10000に固定し、 $N$ を増やしていく。青いヒストグラムは $K$ 個の $\bar{x}$ （平均）の分布を表し、赤い点線はそれぞれの $N$ のときの正規分布

$\begin{align*} p(\bar{x})=\mathcal{N}(\bar{x};\mu,\frac{\sigma^2}{N}) \end{align*}$

である。

図2

$N$ が大きいほど青いヒストグラムの分布は正規分布（赤い点線）に近づいていく様子が分かる。 $N=10$ 程度でほぼ正規分布と一致している。今回使用したソースコードはここにある。

中心極限定理はどこで役立つのか

　ここまでの説明で中心極限定理については理解できたと思う。次にこれが役に立つ場面を紹介する。良く例に使われる日本人の身長について考えよう。日本人の身長を変数 $x$ 、その確率分布を $p(x)$ とおく。 $p(x)$ の真の形状は不明なので（全日本人の身長を測るわけにはいかない）、 $p(x)$ の平均と分散をそれぞれ $\mu$ と $\sigma^2$ という未知数で与えておく。次に $10000(=N)$ 人のサンプリングを行う。

$\begin{align*} D=\{x_1,\cdots,x_N\} \end{align*}$

これらの平均を計算し $\bar{x}$ を求める。このようなサンプリングから得られる平均は、中心極限定理によると次の正規分布

$\begin{align*} p(\bar{x})=\mathcal{N}(\bar{x};\mu,\frac{\sigma^2}{N}) \end{align*}$

に従うのであった。従って、正規分布の場合の区間推定を行うと、 $\bar{x}$ は99.7%の確率で

(1) $\begin{align*} \mu-3\frac{\sigma}{\sqrt{N}}\leq \bar{x} \leq \mu+3\frac{\sigma}{\sqrt{N}} \end{align*}$

を満たすと結論できる。今回は区間推定についての説明は割愛するが、例えば68–95–99.7則を見てほしい。この式を変形すると

$\begin{align*} \bar{x}-3\frac{\sigma}{\sqrt{N}}\leq \mu \leq \bar{x}+3\frac{\sigma}{\sqrt{N}} \end{align*}$

を得る。さらに、統計学の教えるところによると、 $N$ が十分に大きい場合、真の分散 $\sigma^2$ はサンプリングデータから計算できる不偏分散 $s^2$ に置き換えることができる（不偏分散については以前こちらの記事「標本分散と不偏分散」で解説した）。従って

(2) $\begin{align*} \bar{x}-3\frac{s}{\sqrt{N}}\leq \mu \leq \bar{x}+3\frac{s}{\sqrt{N}} \end{align*}$

を得る。上式において未知数は真の平均 $\mu$ だけである。すなわち、真の平均 $\mu$ を99.7%の確率で区間推定できたことになる。このような計算が可能なのは、中心極限定理により $\bar{x}$ が正規分布に従い、その結果、式(1)が成り立つためである。式(2)で $N\rightarrow\infty$ にすると真の平均 $\mu$ が存在する区間が狭くなっていくことが分かる。