この記事はWathematica Advent Calendar 2021の10日目の記事です！
Wathematica Advent Calendar 2021
こんにちは、えぽぱかです。
今回は統計学において基本的で重要な定理である中心極限定理の収束スピードについて考えていきたいと思います！前提知識として数Bの「確率分布と統計的な推測」およびテイラー展開に関する知識を仮定します。
さて、まずは今回扱う中心極限定理のふわっとした主張について見ていきましょう！
定理（ふわふわ中心極限定理）
平均が$μ$, 分散が$σ^2$の確率分布からランダムに$n$個のデータを取ってきたとき, それら$n$個の平均がつくる確率分布は、十分大きい$n$で平均$μ$, 分散$\frac{σ^2}{n}$の正規分布に近似的に従う。
たとえばサイコロを1回振ったときに出た目がつくる確率分布を考えてみましょう。
マトモなサイコロを振ると1から6の目が等確率で出ます。
なので、この確率分布の平均は$(1+2+3+4+5+6)\cdot\frac{1}{6}=3.5$で、分散は$\sum_{i=1}^6(i-3.5)^2 = \frac{35}{12}$となります。
それではこのサイコロを5回振ってみてその出た目の平均を考えてみましょう。この平均をとる操作を数値シミュレーションで何回も繰り返し、求めた平均たちをヒストグラムで表してみると次のようになりました。
実線部分は平均$3.5$, 分散$\frac{35}{12} \cdot \frac{1}{5} = \frac{7}{24}$ の正規分布を表しています。
するとあら不思議、5回の平均たちは正規分布に従っているように見えます！
中心極限定理の主張というのは、サイコロの出た目の確率分布だけではなく、どんな確率分布でも1たくさん平均をとってしまえば正規分布という確率分布に近似できるのだ、ということを言っているのです。
正規分布という特定の確率分布の振る舞いについてちゃんと調べておけばよくわからない確率分布についてなにか問題を考えるときにうまく応用できそうな感じがあってとてもうれしいですね！
でもちょっと待ってください、今の中心極限定理の主張では「十分大きいnだったらn個のデータの平均は正規分布に近似的に従う」と述べました。
サイコロの例ではn=5のとき十分正規分布に従っているように見えましたが、他の確率分布でも5個程度の平均を考えればいい感じに正規分布に近付くと考えてしまってよいのでしょうか？
いいえ、そこまで都合のいいお話ではないのです。しかし、これ以上議論するには中心極限定理をしっかり数学的に見ていく必要があります。
数学的にカッチリした中心極限定理の主張を見るために、前提となる概念について準備しておきましょう。
まず、実数値をとる確率変数$X$について分布関数$F_X(x)$は
$\ \ \ \ \ \ \ \ F_X(x) = P(X\leq x)\ \ \ \ (X$が$x$以下となる確率)2
と定義されます。確率変数が定める確率分布に対して分布関数はただ一つ存在します。
それでは中心極限定理の主張について見ていきましょう。
定理（中心極限定理）
平均が$\mu$, 分散が$\sigma^2$の確率分布からランダムに$n$個のデータ$X_1,\,X_2,\,\ldots$をとってきて, その平均を$\bar{X_n}$とする.
また, 確率変数$X$は平均$0$, 分散$1$の正規分布(標準正規分布)に従うものとする. このとき,
$\displaystyle \ \ \ \ \ \ \ \ \lim_{n\to\infty} P(\frac{\bar{X_n} - \mu}{\frac{\sigma}{\sqrt{n}}} \leq x) = P(X\leq x) = F_X(x)$
が成り立つ.
たくさんデータを持ってきて平均をとれば正規分布に近付くよーって主張が数学的に表現できました。 ここで証明の都合上データの平均$\bar{X_n}$を平均が0, 分散が1の確率分布に乗るように標準化しています。
上の定理が証明できれば$\bar{X_n}$は十分大きいnで平均$μ$, 分散$\frac{σ^2}{n}$の正規分布に近似的に従うことが分かります。
正規分布に近づくスピードについて考察するには中身の証明を見ていく必要があります。証明の前にもう一つ新たな概念の整理が必要です。
確率変数$X$について、$e^{itX}$の期待値である$E[e^{itX}]$を特性関数と呼び、
$\displaystyle \ \ \ \ \ \ \ \ \phi_X(t) := E[e^{itX}]$
と表すことにします。
いきなりよく分からない概念が飛び出してきましたが焦らないでください。実はこの特性関数は確率分布に関してとてもうれしい性質を持っているのです。
1つ目に、分布関数が決まれば特性関数は一意的に決まります！
すなわち、分布関数はあらゆる確率分布に対してただ一つ定まるので、特性関数もあらゆる確率分布に対してただ一つ定まります。
そう、特性関数を調べれば確率分布も調べたことになるのです！3
2つ目に、確率変数の列$X_1,\,X_2,\,\ldots$の特性関数$\phi_{X_n}(t)$を考えてみます。
任意の点$t\in\mathbb{R}$で
$\ \ \ \ \ \ \ \ \displaystyle \lim_{n\to\infty}\phi_{X_n}(t) = \phi_X(t)$
となるとき、
$\phi_{X_n}(t)$に対応する分布関数$F_{X_n}$も$n\to\infty$で$F_X$に収束します！4
これをレヴィの連続性定理と呼びます。5
ふわっと説明すると、nを限りなく大きくしたとき特性関数$\phi_{X_n}(t)$が$\phi_{X}(t)$に限りなく近づいていくなら、
$X_n$の確率分布も$X$の確率分布に近付いていくということですね！
とりあえず、特性関数という概念があって、結構うれしい性質があるよーということを頭に入れておいてください！
それでは、中心極限定理の証明に入っていきましょう。
確率変数$Z$について$E[Z]=\mu,\, V[Z]=\sigma^2$と仮定する.
そして、$Z_1,\, \ldots ,\, Z_n$を$Z$と同じ確率分布からランダムに$n$個とってきたデータとする.
ここで, $X_i = \frac{Z_i - \mu}{\sigma}$として$Z_i$を標準化すると$X_i$は平均0, 分散が1の確率分布に乗る.
また, n個の$Z_i$の平均をそれぞれ$\bar{Z_n}$とおくと,
$\displaystyle \ \ \ \ \ \ \ \ \frac{\bar{Z_n}-\mu}{\frac{\sigma}{\sqrt{n}}} = \frac{X_1 +\cdots + X_n}{\sqrt{n}}$
が簡単な計算により成り立つと分かる. したがって右辺の特性関数が$n\to\infty$で標準正規分布の特性関数に一致することを示せばよい.
$Y_1,\, \ldots ,\, Y_n$を標準正規分布からランダムに$n$個とってきたデータとする。
また,$Z_1,\, \ldots ,\, Z_n$と$Y_1,\, \ldots ,\, Y_n$は独立であるとする. つまり$X_1,\, \ldots ,\, X_n$と$Y_1,\, \ldots ,\, Y_n$も独立となる.
$f(X)=e^{itX}$とおく. すなわち$X$の特性関数$\phi_X(t)$は$E[f(X)]$と表せる.
任意に$t\in\mathbb{R}$をとる.
Taylorの定理より
$\ \ \ \ \ \ \ \ \displaystyle f\left(A+\frac{H}{\sqrt{n}}\right) = f(A) + \frac{1}{\sqrt{n}} f^\prime (A)H + \frac{1}{2n} f^{\prime\prime}(A)H^2 + \mathrm{O}(\frac{1}{n\sqrt{n}})\ \ \ \ (n\to\infty)$
が成り立つ.
$\ \ \ \ \ \ \ \ \displaystyle A_k = \frac{X_1+\cdots +X_{k-1}+Y_{k+1}+\cdots +Y_n}{\sqrt{n}}$
とおくと, $A_k,\, X_k,\, Y_k$は独立である.
$\ \ \ \ \ \ \ \ \displaystyle f\left(\frac{X_1+\cdots +X_k + Y_{k+1}+\cdots +Y_n}{\sqrt{n}}\right) - f\left(\frac{X_1+\cdots +X_{k-1} + Y_k +\cdots +Y_n}{\sqrt{n}}\right) = f\left(A_k + \frac{X_k}{\sqrt{n}} \right) - f\left(A_k + \frac{Y_k}{\sqrt{n}} \right)$
$\displaystyle \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \frac{1}{\sqrt{n}}f^\prime (A_k)(X_k-Y_k) + \frac{1}{2n}f^{\prime\prime}(A_k)(X_k^2 - Y_k^2) + \mathrm{O}(\frac{1}{n\sqrt{n}}).\ \ \ \ (n\to\infty)$
ここで, $A_k,\, X_k,\, Y_k$は独立なので$E[f^{\prime}(A_k)(X_k-Y_k)]$ $=$ $E[f^{\prime}(A_k)]$ $(E[X_k]-E[Y_k])$ $=$ $0$,
$E[f^{\prime\prime}(A_k)(X_k^2-Y_k^2)]=0.$
よって, 上式の両辺の期待値をとって
$\ \ \ \ \ \ \ \ \displaystyle E\left[f\left(\frac{X_1+\cdots +X_k + Y_{k+1}+\cdots +Y_n}{\sqrt{n}}\right)\right] - E\left[f\left(\frac{X_1+\cdots +X_{k-1} + Y_k +\cdots +Y_n}{\sqrt{n}}\right)\right] = \mathrm{O}(\frac{1}{n\sqrt{n}}).\ \ (n\to\infty)$
この式を$k=1,\,\ldots,\, n$について足し上げると,$\displaystyle f\left(A_k + \frac{Y_k}{\sqrt{n}}\right) = f\left(A_{k-1} + \frac{X_{k-1}}{\sqrt{n}}\right)$であるから, 打ち消し合って$k=n$のときの第1項と$k=1$のときの第2項のみが残って
$\ \ \ \ \ \ \ \ \displaystyle E\left[f\left(\frac{X_1 + \cdots + X_n}{\sqrt{n}}\right)\right] - E\left[f\left(\frac{Y_1 + \cdots + Y_n}{\sqrt{n}}\right)\right] = \mathrm{O}(\frac{1}{\sqrt{n}})\to 0.\ \ \ \ (n\to\infty)$
$Y_i$は平均$0$, 分散$1$の標準正規分布に従うので, 正規分布の再生性と簡単な計算より$\displaystyle \frac{Y_1+\cdots +Y_n}{\sqrt{n}}$も標準正規分布に従うから,
標準正規分布の特性関数を$\phi(t)$とおくと
$\displaystyle \ \ \ \ \ \ \ \ E\left[f\left(\frac{Y_1 + \cdots + Y_n}{\sqrt{n}}\right)\right] = \phi(t)$
がnに依存せず成立する. したがって
$\displaystyle \ \ \ \ \ \ \ \ \lim_{n\to\infty} E\left[f\left(\frac{\bar{Z_n}-\mu}{\frac{\sigma}{\sqrt{n}}}\right)\right] = \lim_{n\to\infty} E\left[f\left(\frac{X_1 + \cdots +X_n}{\sqrt{n}}\right)\right] = \phi(t).$
よってレヴィの連続性定理より$\frac{\bar{Z_n}-\mu}{\frac{\sigma}{\sqrt{n}}}$の分布関数は$n\to\infty$で標準正規分布の分布関数に一致する.
証明ができました。収束スピードについて調べるために$\mathrm{O}(\frac{1}{n\sqrt{n}})$の部分を計算してみましょう。$E[X^3]$が有限であると仮定します。
$\displaystyle \ \ \ \ \ \ \ \ \mathrm{O}(\frac{1}{n\sqrt{n}}) = \frac{E[f^{\prime\prime\prime}(A_k)]}{6n\sqrt{n}}(E[X_k^3] - E[Y_k^3]) + \mathrm{O}(\frac{1}{n^2})$
$\displaystyle \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \frac{E[f^{\prime\prime\prime}(A_k)]}{6n\sqrt{n}}E[\left(\frac{Z_k - \mu}{\sigma}\right)^3] + \mathrm{O}(\frac{1}{n^2}). \ \ \ \ (\because Y_k$は標準正規分布に従うので$E[Y_k^3]=0)$
$\displaystyle \ \ \ \ \ \ \ \ \left\lvert \frac{E[f^{\prime\prime\prime}(A_k)]}{6n\sqrt{n}}E[\left(\frac{Z_k - \mu}{\sigma}\right)^3]\right\rvert = \frac{\lvert t\rvert^3}{6n\sqrt{n}}\left\lvert \frac{E[(Z_k - \mu)^3]}{\sigma^3}\right\rvert.$
よって収束する速さは$\left\lvert \frac{E[(Z_k - \mu)^3]}{\sigma^3}\right\rvert$の大きさでおおよそ決まることが分かります.
$\frac{E[(Z_k - \mu)^3]}{\sigma^3}$という量には実は名前がついており、歪度と呼びます。歪度は分布の歪み具合を端的に表す量だと考えることができ、歪度0の確率分布は左右対称であり、歪度の絶対値が大きい確率分布はあまり対称性のない分布であると言えます。
すなわち、大まかに言うと中心極限定理の収束スピードは元の確率分布の歪度によって決まり、左右対称な確率分布ほど収束スピードが速いのです！ 序盤で示した一様分布は左右対称な分布なのでn=5程度ですぐに標本平均の分布が正規分布に近くなったということなのですね。
それでは最後に、いろんな確率分布を数値シミュレーションしてみて中心極限定理の収束の具合を見ていきましょう！
シミュレーションのコードはMathematicaで書きました。大学にもよりますが学生は無料でインストールして使うことができます！使ったコードはこちらです。
n = 5; (* 何個の平均を考えるか *)
m = 10^8; (* 合計で何個のデータをとってくるか *)
dist =
DiscreteUniformDistribution[{1, 6}]; (* 母集団の確率分布 *)
mu =
Mean[dist]; v = Variance[dist];
data = RandomVariate[dist, m];
smdata = Mean[Partition[data, m/n]];
hist = Histogram[smdata, Automatic, "ProbabilityDensity"]; (* ヒストグラムがうまく表示されないときはAutomaticの部分を100とかに書き換える*)
pl = Plot[
PDF[NormalDistribution[mu, Sqrt[v/n]], x], {x, mu - 3* Sqrt[v/n],
mu + 3* Sqrt[v/n]}];
Show[pl, hist]
まずはベータ分布$Beta(0.5,\, 0.5)$で試してみましょう。
$Beta(0.5,\, 0.5)$の歪度は0となるのでデータの個数が少なくてもその平均は正規分布に従いそうです。実際にn=5としてシミュレーションを実行してみましょう。
なんか惜しい感じがします。n=10で再度シミュレーションしてみましょう。
いい感じに正規分布に従ってそうですね！
次に指数分布$Ex(1)$を考えてみましょう。
指数分布の歪度は２となるのでデータの個数が少ないとあまり正規分布に近付かなそうです。n=10でシミュレーションしてみましょう。
たしかになんか渋いですね。n=100だとどうでしょうか。
それっぽさは感じますがもう一声ほしくなります。n=500で行ってみましょう。
まあよいでしょう。このように歪度によって中心極限定理の収束スピードは全然変わってきます。
次にラプラス分布$Laplace(0,1)$を考えてみましょう。
ラプラス分布の歪度は０となるのでデータの個数が少なくても正規分布に近付きそうです。さあどうでしょうか。n=10でシミュレーションしてみましょう。
見た感じ良さそうですがよく見ると山の頂上の部分が若干尖っている感じがします。
先ほど中心極限定理の収束する速さはおおよそ歪度で決まると言いましたが、歪度が0のとき同じような議論で収束スピードは$\frac{E[(Z_k - \mu)^4]}{\sigma^4}$で決まることが分かります。
この量を尖度と言いますが、ラプラス分布の尖度は3なのでちょっぴり収束スピードが遅くなるのですね。
しかし、歪度よりは速度に影響してこないのでn=50で大体正規分布に従いそうなことが以下のシミュレーション結果で分かります。
最後にフレシェ分布$Frechet(2.5,1)$を考えてみましょう。
フレシェ分布は形状パラメータが2.5のとき歪度が無限大となります。（期待値と分散は有限値になります)
歪度が無限大となるので先程の収束スピードの議論はうまく適用できませんが中心極限定理は成り立ちます。
直感的には非常に大きいnでないと正規分布に従わなさそうですね。実際n=100でシミュレーションしてみましょう。
全然従ってる感じがしないですね。n=1000ではどうでしょう。
これでもアレですね。n=10000で行っときますか。
まあ耐えてるのかこれは……？とにかく収束するスピードはめちゃくちゃ遅いことが分かります。
いかがでしたか？これにて今回の記事を締めたいと思います！ありがとうございました！
参考文献
- 黒木玄, “確率論入門,” 1 9 2017. [オンライン]. Available: https://genkuroki.github.io/documents/IntroProbability.pdf.
- 鈴木将史, “確率分布の収束とキュムラント,” 1991. [オンライン]. Available: https://core.ac.uk/download/pdf/147576366.pdf.
- 秋田智之, “漸近展開へのいざない,” 24 3 2009. [オンライン]. Available: https://home.hiroshima-u.ac.jp/akita/stat/Edgeworth.pdf.
- 久保川達也, 現代数理統計学の基礎, 共立出版, 2017.