✏️
中心極限定理において, 標準正規分布に早く近付くのはどんな分布か？

2025/05/10に公開
 これは何かタイトルの通り
数学的な厳密さは気にせずに記述する

 表記確率変数Xの特性関数: \varphi_{X}(t) = E \left[ \exp\left( itX \right) \right]

\varphi_{X}(t)のk階微分: \varphi_{X}^{(k)} (t)


 問題設定
X_{k} (k=1,\ldots,n)を, 確率密度関数f_{X}(x)に独立に従う確率変数とする
平均: \mu = E\left[ X \right]
分散: \sigma^{2} = E\left[ \left( X - \mu \right)^{2} \right]
歪度: \gamma = \frac{E\left[ \left( X - \mu \right)^{3} \right]}{\sigma^{3}}
標本平均: \bar{X} = \frac{1}{n} \sum_{k=1}^{n} X_{k}


Z = \frac{\sqrt{n} \left( \bar{X} - \mu \right)}{\sigma} とおくと, 中心極限定理より n \rightarrow \inftyでZは標準正規分布に分布収束することがわかる
中心極限定理において, 標準正規分布に早く近付くのはどんな分布か？について考えたい

 結論中心極限定理において, 標準正規分布に早く近付くのはどんな分布か？
先に答えだけ述べると, 歪度が0に近い 分布である.

 数式で確認中心極限定理を証明する中で, 正規分布への収束の早さが歪度に依存することが確認できる.
Zの特性関数\varphi(t)は

\begin{aligned}
\varphi(t)
&= E \left[ \exp \left( itZ \right) \right] = E \left[ \exp \left( it \frac{\sqrt{n} \left( \bar{X} - \mu \right)}{\sigma} \right) \right] \\
&= E \left[ \exp \left( i \frac{t}{\sqrt{n} \sigma} \sum_{k=1}^{n} X_{k} \right) \right] \exp \left( - \frac{i \sqrt{n} \mu t}{\sigma} \right) \\
&= \left\{ \varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) \right\}^{n} \exp \left( - \frac{i \sqrt{n} \mu t}{\sigma} \right)
\end{aligned}
となるが, 前半の\left\{ \varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) \right\}^{n}の項をTaylor展開を用いて評価していくことで

\begin{aligned}
\varphi(t)
&= \left\{ \varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) \right\}^{n} \exp \left( - \frac{i \sqrt{n} \mu t}{\sigma} \right) \\
&= \exp \left( n^{\frac{1}{2}} \frac{i \mu t}{\sigma} - \frac{t^{2}}{2} - n^{-\frac{1}{2}} \frac{i \gamma t^{3}}{6} + \cdots \right) \exp \left( - \frac{i \sqrt{n} \mu t}{\sigma} \right) \\
&= \exp \left( - \frac{t^{2}}{2} - n^{-\frac{1}{2}} \frac{i \gamma t^{3}}{6} + \cdots \right)
\end{aligned}
という式が得られる（計算の詳細は後述）.
この式から, Zの特性関数はn \rightarrow \inftyで標準正規分布の特性関数\exp \left( - \frac{t^{2}}{2} \right)に収束することがわかり, その際の近付き方は歪度 \gamma が0に近いほど早いことがわかる.
計算の詳細\varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) のTaylor展開を考えると

\begin{aligned}
& \varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) \\
=& \varphi_{X} \left( 0 \right) + \varphi_{X}^{(1)} \left( 0 \right) \left( \frac{t}{\sqrt{n} \sigma} \right) + \frac{1}{2!} \varphi_{X}^{(2)} \left( 0 \right) \left( \frac{t}{\sqrt{n} \sigma} \right)^{2} + \frac{1}{3!} \varphi_{X}^{(3)} \left( 0 \right) \left( \frac{t}{\sqrt{n} \sigma} \right)^{3} + \cdots \\
=& 1 + E \left[ iX \right] \left( \frac{t}{\sqrt{n} \sigma} \right) + \frac{1}{2!} E \left[ \left( iX \right)^{2} \right] \left( \frac{t}{\sqrt{n} \sigma} \right)^{2} + \frac{1}{3!} E \left[ \left( iX \right)^{3} \right] \left( \frac{t}{\sqrt{n} \sigma} \right)^{3} + \cdots \\
=& 1 + \frac{i \mu t}{\sqrt{n} \sigma} - \frac{E \left[ X^{2} \right] t^{2}}{2 n \sigma^{2}} - \frac{i E \left[ X^{3} \right] t^{3}}{6 n^{\frac{3}{2}} \sigma^{3}}  + \cdots
\end{aligned}
となるため, このn乗は

\begin{aligned}
& \left\{ \varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) \right\}^{n} \\
=& \exp \left[ n \log \varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) \right] \\
=& \exp \left[ n \log \left\{ 1 + \frac{i \mu t}{\sqrt{n} \sigma} - \frac{E \left[ X^{2} \right] t^{2}}{2 n \sigma^{2}} - \frac{i E \left[ X^{3} \right] t^{3}}{6 n^{\frac{3}{2}} \sigma^{3}}  + \cdots \right\} \right]
\end{aligned}
となる. \logの項を \log (1 + x) = x - \frac{1}{2} x^{2} + \frac{1}{3} x^{3} - \cdotsのTaylor展開で展開すると

\begin{aligned}
& \log \left\{ 1 + \frac{i \mu t}{\sqrt{n} \sigma} - \frac{E \left[ X^{2} \right] t^{2}}{2 n \sigma^{2}} - \frac{i E \left[ X^{3} \right] t^{3}}{6 n^{\frac{3}{2}} \sigma^{3}}  + \cdots \right\} \\
=& \left\{ \frac{i \mu t}{\sqrt{n} \sigma} - \frac{E \left[ X^{2} \right] t^{2}}{2 n \sigma^{2}} - \frac{i E \left[ X^{3} \right] t^{3}}{6 n^{\frac{3}{2}} \sigma^{3}}  + \cdots \right\} \\
 & - \frac{1}{2} \left\{ \frac{i \mu t}{\sqrt{n} \sigma} - \frac{E \left[ X^{2} \right] t^{2}}{2 n \sigma^{2}} - \frac{i E \left[ X^{3} \right] t^{3}}{6 n^{\frac{3}{2}} \sigma^{3}}  + \cdots \right\}^{2} \\
 & + \frac{1}{3} \left\{ \frac{i \mu t}{\sqrt{n} \sigma} - \frac{E \left[ X^{2} \right] t^{2}}{2 n \sigma^{2}} - \frac{i E \left[ X^{3} \right] t^{3}}{6 n^{\frac{3}{2}} \sigma^{3}}  + \cdots \right\}^{3} \\
 & + \cdots \\
=& n^{-\frac{1}{2}} \frac{i \mu t}{\sigma} - n^{-1} \frac{ \left( E \left[ X^{2} \right] - \mu^{2} \right) t^{2}}{2 \sigma^{2}} - n^{-\frac{3}{2}} \frac{i \left( E \left[ X^{3} \right] - 3 E \left[ X^{2} \right] \mu + 2 \mu^{3} \right) t^{3}}{6 \sigma^{3}} + \cdots \\
=& n^{-\frac{1}{2}} \frac{i \mu t}{\sigma} - n^{-1} \frac{\sigma^{2} t^{2}}{2 \sigma^{2}} - n^{-\frac{3}{2}} \frac{i E \left[ \left( X - \mu \right)^{3} \right] t^{3}}{6 \sigma^{3}} + \cdots \\
=& n^{-\frac{1}{2}} \frac{i \mu t}{\sigma} - n^{-1} \frac{t^{2}}{2} - n^{-\frac{3}{2}} \frac{i \gamma t^{3}}{6} + \cdots \\
\end{aligned}
となるため,

\begin{aligned}
& \left\{ \varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) \right\}^{n} \\
=& \exp \left[ n \log \left\{ 1 + \frac{i \mu t}{\sqrt{n} \sigma} - \frac{E \left[ X^{2} \right] t^{2}}{2 n \sigma^{2}} - \frac{i E \left[ X^{3} \right] t^{3}}{6 n^{\frac{3}{2}} \sigma^{3}}  + \cdots \right\} \right] \\
=& \exp \left[ n \left\{ n^{-\frac{1}{2}} \frac{i \mu t}{\sigma} - n^{-1} \frac{t^{2}}{2} - n^{-\frac{3}{2}} \frac{i \gamma t^{3}}{6} + \cdots \right\} \right] \\
=& \exp \left( n^{\frac{1}{2}} \frac{i \mu t}{\sigma} - \frac{t^{2}}{2} - n^{-\frac{1}{2}} \frac{i \gamma t^{3}}{6} + \cdots \right)
\end{aligned}
となる.
以上より, Zの特性関数\varphi(t)は

\begin{aligned}
\varphi(t)
&= \left\{ \varphi_{X} \left( \frac{t}{\sqrt{n} \sigma} \right) \right\}^{n} \exp \left( - \frac{i \sqrt{n} \mu t}{\sigma} \right) \\
&= \exp \left( n^{\frac{1}{2}} \frac{i \mu t}{\sigma} - \frac{t^{2}}{2} - n^{-\frac{1}{2}} \frac{i \gamma t^{3}}{6} + \cdots \right) \exp \left( - \frac{i \sqrt{n} \mu t}{\sigma} \right) \\
&= \exp \left( - \frac{t^{2}}{2} - n^{-\frac{1}{2}} \frac{i \gamma t^{3}}{6} + \cdots \right)
\end{aligned}
となる.
（計算の詳細終わり）

 実際の分布で確認Zの分布を書き下すのが容易な分布で, 本当に（元の分布の）歪度が小さいほどZの標準正規分布への近付き方が早いのかを確認してみる.

 ベルヌーイ分布ベルヌーイ分布Bernoulli\left( p \right)の基本情報
確率関数: P \left( X = x \right) = p^{x} \left( 1 - p \right)^{1-x}
平均: \mu = p
分散: \sigma^{2} = p (1 - p)
歪度: \gamma = \frac{1 - 2p}{\sqrt{p (1 - p)}}


Z = \frac{\sqrt{n} \left( \bar{X} - \mu \right)}{\sigma}の分布関数
基本的には, 二項分布の分布関数と同じ形
分布関数を書き下すと以下のようになる:


F_{Z} \left( z \right) = \sum_{k=0}^{\lfloor \sqrt{n p (1 - p)} z + np \rfloor} \binom{n}{k} p^{k} \left( 1 - p \right)^{n - k}
ただし, \lfloor s \rfloorはsを越えない最大の整数とする. また, \lfloor \sqrt{n p (1 - p)} z + np \rfloorの値が0未満の場合は分布関数の値は0, nより大きい場合は分布関数の値は1となる.
Zの分布関数の導出まず, Zは

\begin{aligned}
Z
&= \frac{\sqrt{n} \left( \bar{X} - \mu \right)}{\sigma} 
= \frac{\sqrt{n} \left( \bar{X} - p \right)}{\sqrt{p (1 - p)}} \\
&= \frac{\sum_{k=1}^{n} X_{k} - n p}{\sqrt{n p (1 - p)}}
\end{aligned}
となる.
ここで, ベルヌーイ分布の変数の和S = \sum_{k=1}^{n} X_{k} は二項分布Binom \left( n, p \right)に従うため, このSの分布関数は

F_{S} \left( s \right) = P \left( S \le s \right) = \sum_{k=0}^{\lfloor s \rfloor} \binom{n}{k} p^{k} \left( 1 - p \right)^{n - k}
と書ける. ただし, \lfloor s \rfloorはsを越えない最大の整数とする.
するとZの分布関数は

\begin{aligned}
F_{Z} \left( z \right)
=& = P \left( Z \le z \right) = P \left( \frac{S - n p}{\sqrt{n p (1 - p)}} \le z \right) \\
=& P \left( S \le \sqrt{n p (1 - p)} z + np \right) \\
=& \sum_{k=0}^{\lfloor \sqrt{n p (1 - p)} z + np \rfloor} \binom{n}{k} p^{k} \left( 1 - p \right)^{n - k}
\end{aligned}
となる. ただし, \lfloor \sqrt{n p (1 - p)} z + np \rfloorの値が0未満の場合は分布関数の値は0, nより大きい場合は分布関数の値は1となる.
（Zの分布関数の導出終わり）
サンプルサイズnと歪度\frac{1 - 2p}{\sqrt{p (1 - p)}}を動かし, Zの確率密度関数の形がどうなるのかを観察する. これまでの議論を踏まえると, 歪度\frac{1 - 2p}{\sqrt{p (1 - p)}}が0に近いほど, つまりpが\frac{1}{2}に近いほどZは標準正規分布に早く収束することが期待される.
可視化コードclt_bernoulli.py
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm, binom


def plot_bernoulli():
    ns = [5, 10, 50]
    ps = [0.05, 0.1, 0.5, 0.9, 0.95]
    zs = np.linspace(-3, 3, 101)
    _, axes = plt.subplots(1, len(ns), figsize=(8, 12))
    for fig_index, n in enumerate(ns):
        for p in ps:
            cdfs = binom.cdf(np.floor(zs * np.sqrt(n * p * (1 - p)) + n * p), n, p)
            axes[fig_index].plot(zs, cdfs, label=f"Z(p={p})")
        axes[fig_index].plot(zs, norm.cdf(zs, loc=0, scale=1), label="Normal(0, 1)", color="black")
        axes[fig_index].legend(loc="upper left")
        axes[fig_index].set_title(f"Cumulative Distributions(n={n})")
    plt.tight_layout()
    plt.show()


def main():
    plot_bernoulli()


if __name__ == "__main__":
    main()
nとpを動かし, Zの分布関数と標準正規分布の分布関数を重ねてプロットしてみると, 以下のようになる:
グラフを見る限り, pが\frac{1}{2}に近いほど, つまり歪度\frac{1 - 2p}{\sqrt{p (1 - p)}}が0に近いほどZは標準正規分布に早く収束することが確認できる.

 ガンマ分布ガンマ分布Gamma\left( \alpha, \beta \right)の基本情報
確率密度関数: f(x) = \frac{\beta^{\alpha}}{\Gamma \left( \alpha \right)} x^{\alpha - 1} \exp \left( - \beta x \right)
平均: \mu = \frac{\alpha}{\beta}
分散: \sigma^{2} = \frac{\alpha}{\beta^{2}}
歪度: \gamma = \frac{2}{\sqrt{\alpha}}


Z = \frac{\sqrt{n} \left( \bar{X} - \mu \right)}{\sigma}の確率密度関数
ガンマ分布Gamma\left( n \alpha, \sqrt{n \alpha} \right)を負の方向に\sqrt{n \alpha}だけずらしたもの（導出は後述）
確率密度関数を書き下すと, 以下の通り:


f_{Z}(z) = \frac{\left( n \alpha \right)^{\frac{\alpha}{2}}}{\Gamma \left( n \alpha \right)} \left( z + \sqrt{n \alpha} \right)^{n \alpha - 1} \exp \left\{ - \sqrt{ n \alpha } \left( z + \sqrt{n \alpha} \right) \right\}
Zの確率密度関数の導出まず, Zは

\begin{aligned}
Z
&= \frac{\sqrt{n} \left( \bar{X} - \mu \right)}{\sigma} 
= \frac{\sqrt{n} \left( \bar{X} - \frac{\alpha}{\beta} \right)}{\sqrt{\frac{\alpha}{\beta^{2}}}} \\
&= \frac{\beta}{\sqrt{n \alpha}} \sum_{k=1}^{n} X_{k} - \sqrt{n \alpha} \\
\end{aligned}
となる. 第2項は定数なので, Zの分布は「第1項の分布を負の方向に\sqrt{n \alpha}だけ平行移動した分布」となる.
そこで, 第1項がどんな分布に従うかを考える. ガンマ分布の再生性より, \sum_{k=1}^{n} X_{k}はガンマ分布Gamma\left( n \alpha, \beta \right) に従う. 第1項はそれを\frac{\beta}{\sqrt{n \alpha}}倍したものであるため, ガンマ分布Gamma\left( n \alpha, \beta \times \frac{\sqrt{n \alpha}}{\beta} \right) = Gamma\left( n \alpha, \sqrt{n \alpha} \right)に従う.
以上のことから, Zの従う分布の確率密度関数は

f_{Z}(z) = \frac{\left( n \alpha \right)^{\frac{\alpha}{2}}}{\Gamma \left( n \alpha \right)} \left( z + \sqrt{n \alpha} \right)^{n \alpha - 1} \exp \left\{ - \sqrt{ n \alpha } \left( z + \sqrt{n \alpha} \right) \right\}
となる.
（Zの確率密度関数の導出終わり）
サンプルサイズnと歪度\frac{2}{\sqrt{\alpha}}を動かし, Zの確率密度関数の形がどうなるのかを観察する. これまでの議論を踏まえると, 歪度\frac{2}{\sqrt{\alpha}}が0に近いほど, つまり\alphaが大きいほどZは標準正規分布に早く収束することが期待される. ただし, 確率密度関数の式からわかるように, 以下の点に留意する必要がある:

nと\alphaどちらを変化させても, Zの確率密度関数の形への影響の仕方は変わらない

Zの確率密度関数は\betaに依存しない
可視化コードclt_gamma.py
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm, gamma


def plot_gamma():
    alphas = [1, 2, 4]
    ns = [5, 10, 50]
    zs = np.linspace(-3, 3, 1001)
    _, axes = plt.subplots(1, len(ns), figsize=(8, 12))
    for fig_index, n in enumerate(ns):
        for alpha in alphas:
            gamma_ps = gamma.pdf(zs, a=n * alpha, loc=-np.sqrt(n * alpha), scale=1/np.sqrt(n * alpha))
            axes[fig_index].plot(zs, gamma_ps, label=f"Z(alpha={alpha})")
        axes[fig_index].plot(zs, norm.pdf(zs, loc=0, scale=1), label=f"Normal(0, 1)", color="black")
        axes[fig_index].set_title(f"Probability Densities (n={n})")
        axes[fig_index].legend(loc="upper right")
    plt.tight_layout()
    plt.show()


def main():
    plot_gamma()


if __name__ == "__main__":
    main()
nと\alphaを動かし, Zの確率密度関数と標準正規分布の確率密度関数を重ねてプロットしてみると, 以下のようになる:
グラフを見る限り, \alphaが大きいほど, つまり歪度\frac{2}{\sqrt{\alpha}}が0に近いほどZは標準正規分布に早く収束することが確認できる.

 所感個人的には, 中心極限定理の証明を追いかける中で歪度が自然に登場してきたことに感動した
歪度については, 平均, 分散と比べると影が薄い印象がある
統計学の入門書だと, 中心極限定理に関しては特性関数が\exp \left( - \frac{1}{2} t^{2} \right)に収束することを示して終わり, という印象が強いが, n^{- \frac{1}{2}}の項まで見ることで歪度が登場してくるのはなかなか興味深いなと感じた
これは何か

表記

問題設定

結論

数式で確認

実際の分布で確認

ベルヌーイ分布

ガンマ分布

所感

Discussion
