カイ二乗分布、t分布、F分布の導出


#概要

以下のノートでは、$\chi^2$分布、t分布、F分布を導出したいと思います。これらの分布は、中級統計学の授業で出てきますが、導出はあまり紹介されないことが多く、統計学を学ぶ上での（あるいは方便として統計検定を学ばなければいけない時の）障害になっているのではないかと思います。

これらの分布は、多少面倒ではありますが、比較的見通しの良い方法で導かれることを主張したいと思います。母集団の分布が正規分布であると仮定し、観測値$X_i$が得られたとき、ベイズの定理を用いると、パラメーターの確率分布$P(\mu, \sigma | X_i )$が得られます。ここで、$\mu$は平均、$\sigma$は標準偏差です。このとき、$\chi^2$分布、t分布、F分布は、この分布を用いて求める変数を固定し、その他のすべての変数について積分することにより得られます。$\chi^2$分布は変数を$\sigma^2$のみ残した時の分布であり、t分布は$\mu$の分布、F分布は$\sigma_y/\sigma_y$（2つの正規分布の標準偏差の比）の分布となります。

ただし、自分は統計学の専門ではないため、以下の導出に問題や疑問等がありましたら、（特にご専門の方から）コメント頂けると有難いです。

#準備

ベイズの定理を利用するためには，事前分布を求める必要があります。我々の目的の分布を導くためにはどうしても事前分布が重要になりますが，我々はいわゆる無情報事前分布と呼ばれるものを用いる事にします。（ジェフリーの事前分布と呼ばれるものです）
以下の準備は一般的でない内容を含みますが、ジェフリーの事前分布という、なるべく事前分布が情報を持たないような事前分布を採用して導出を行うことだけ認めて頂ければ大丈夫です。

しかし，事前分布の「無情報性」を判断するには，どのような基準があるのでしょうか？
それは、元のパラメーター$\mu$、$\sigma$を変数変換することにより得られた変数組に関して、
その空間が「均一」であることを課すことによって得られます。

ここでは、$\sigma$を変数変換し、$f(\sigma)$という新たな変数を用いるとしましょう。このとき、$f$を一定区間ごとに区切ったときの確率分布の変化が、一定の変化率で起こるようなパラメーター$f$がとれれば、その空間は均一にパラメトライズできたと思えるのではないでしょうか。
実は、確率分布の変化の大きさはカルバック・ライブラ距離という物で与えられますが、$f$を$f+df$に変化させたときの確率分布の変化の大きさが常に一定であれば、新しいパラメーター$f$について確率分布の変化は一定だと思えるので、そのような$f$について$df=\rho d\sigma$となる$\rho$が事前分布としてふさわしいと考える事ができます。

よって、ここでは、$df$の変化に対するカルバック・ライブラ距離が常に一定になるという条件より、（分布は正規分布であると仮定して、$X$, $\mu$, $\sigma$を用いて$P(X, \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^2}}\exp\left( -\frac{(X-\mu)^2}{2 \sigma^2} \right)$と書けるものとします）

```math
\begin{equation}
\int P(f) \  \log \frac{P(f+df)}{P(f)} dX = C
\end{equation}
```
という条件が得られます。
左辺を展開し、

```math
\begin{equation}
\begin{split}
\int P(f) \log \frac{P(f+df)}{P(f)} dX &
\sim \int P(f) \left\{ df \partial_f (\log P(f))+ \frac{1}{2} df^2 \partial^2_f(\log P(f)) \right\}\\
& \sim df^2 \left( \frac{d\sigma}{df} \right)^2 \frac{1}{\sigma^2}
\end{split}
\end{equation}
```
1行目から2行目の式変形では、愚直に計算を行い、dfの1次の項は消え、2次の項は上式にあるように(今は$df$を等間隔と考えているので)$(d\sigma / df) (1/\sigma^2)$に比例する項が得られます。
よってこの項が一定であることから$df \sim d\sigma / \sigma$が得られ、事前分布$\rho$を $\rho = 1/\sigma$ととるのが自然であることがわかります。

#カイ二乗分布の導出
$\chi^2$分布は、母集団の分布が正規分布である場合に、標準偏差$\sigma$を変数変換した変数$Y$に関する分布となります。ここで、標準偏差は$Y\equiv \frac{n-1}{\sigma^2}s$ととります。ここで、$\sigma$：母集団分布の未知の標準偏差、n：観測されたサンプル数、$s^2$：観測された分散であり、$s^2 = \frac{1}{n-1} \sum (X_i - \bar{X})^2$、ただし$\bar{X}\equiv\frac{1}{n}\sum X_i$であるとします。
ここで、$\sigma$を$Y$と変数変換したのは恣意的にも見えますが、これは、のちに得られる式が既知の積分の形（ガンマ関数）になるように調整しただけだということが見て頂けると思います。
ベイズの定理より、パラメーター$\mu$、$\sigma$に関する分布は

```math
\begin{equation}
P(\mu, \sigma | X_i) \propto P(X_i|\mu,\sigma) P(\mu)P(\sigma)    
\end{equation}
```
と書けます。
定数倍の違いは気にしないことにすると、$\mu, \sigma$に関する積分は次のように書けます。

```math
\begin{equation}
\mathcal{Z} = \int P(\mu, \sigma | X_i) d\mu d\sigma \propto \int P(X_i|\mu,\sigma) \frac{1}{\sigma} d\mu d\sigma.
\end{equation}
```

母集団の分布は正規分布であると仮定しているため、
$P(X_i|\mu, \sigma)= \prod_i \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left( -\frac{(X_i-\mu)^2}{2\sigma^2} \right)$
であり、さらに式変形して$\mu$について積分すると、

```math
\begin{equation}
\begin{split}
\mathcal{Z} & \sim \int d\mu d\sigma (2\pi)^{-\frac{n}{2}}\sigma^{-n} \exp \left(  -\frac{1}{2\sigma^2}(n\mu^2-2n\bar{X}\mu+\sum_i X_i^2) \right) \frac{1}{\sigma}\\
            & \sim \int d\mu d\sigma (2\pi)^{-\frac{n}{2}}\sigma^{-n} \exp \left( -\frac{n}{2\sigma^2}(\mu-\bar{X})^2+\frac{n}{2\sigma^2}(\bar{X}^2-\frac{1}{n}\sum_i X_i^2) \right) \frac{1}{\sigma}\\
            & \sim \int d\sigma (2\pi)^{-\frac{n}{2}+1}n^{-\frac{1}{2}}\sigma^{-n+1} \exp\left( -\frac{(n-1)}{2\sigma^2}s^2 \right)\frac{1}{\sigma}
\end{split}
\end{equation}
```

となります。ここで2行目から3行目の式変形では不偏標準偏差の表式より、$s^2 = \frac{n}{n-1}(\frac{1}{n}\sum X_i^2-\bar{X}^2)$を用いました。

さらに、この積分をガンマ関数の式に帰着させるため、新たな変数を$Y=\frac{n-1}{\sigma^2}s^2$ととると、その微分要素は$dY= -2\frac{Y}{\sigma}d\sigma$となり、全ての定数倍の表式を無視して、さらに自由度の表式$k \equiv n-1$を用いると、最終的に

```math
\begin{equation}
\mathcal{Z} \sim \int dY Y^{\frac{k}{2}-1}\exp\left( -\frac{Y}{2} \right). 
\end{equation}
```

が得られます。
これは我々が求めようとしていた$\chi^2$分布の式です。

#t分布の導出
$\chi^2$分布の導出と同様に、我々はここでは、与えられた観測値が$X_1 \dots X_n$であり、母集団の分布は正規分布であるとして、真の平均の値$\mu$の分布を求めます。
$P(\mu,\sigma|X_i)$についてベイズの定理を用いて、さらに$\sigma$について積分すると、

```math
\begin{equation}
    \begin{split}
        \mathcal{Z}=& \int d\mu d\sigma P(\mu,\sigma|X_i)\\
        \sim & \int d\mu d\sigma P(X_i|\mu,\sigma)P(\mu)P(\sigma)\\
        = & \int d\mu d\sigma \prod_i \left[\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(X_i-\mu)^2}{2\sigma^2} \right) \right]\frac{1}{\sigma}\\
    \end{split}
\end{equation}
```

となります。
1行目の表式と2行目の表式は（規格化のための）定数倍だけ異なります(この定数は$\mu$、$\sigma$に依存します)。$\sigma$に関する積分を実行するため、さらに変数変換を行い、$\frac{1}{2\sigma^2}\equiv\sigma'$とし、このとき微分要素は$d\sigma/\sigma = -d\sigma' /(2\sigma')$となります。すると

```math
\begin{equation}
    \begin{split}
    \mathcal{Z} \sim & \int d\mu d\sigma' \left( -\frac{1}{2\sigma'}\right) (2\pi)^{-\frac{n}{2}}(2\sigma')^{\frac{n}{2}}\exp\left( -\sum_i (X_i-\mu)^2 \sigma' \right)\\
        \sim & \int d\mu d\sigma' \sigma'^{\frac{n}{2}-1} \exp\left( -\sum_i (X_i-\mu)^2 \sigma' \right)
    \end{split}
\end{equation}
```

となります。ここでまた定数倍の違いは無視した表式となっています。さらに$\sigma'$に関して積分すると、

```math
\begin{equation}
    \begin{split}
        \mathcal{Z} \sim \int d\mu ( \sum_i (X_i-\mu)^2 )^{-\frac{n}{2}}.
    \end{split}
\end{equation}
```

となります。
ここでさらに$\sum_i (X_i-\mu)^2=n\mu^2-2\mu n\bar{X} +\sum_i X_i^2=n(\mu-\bar{X})^2 + (n-1)s^2$（ただし $s^2 = \frac{1}{n-1}\sum_i (X_i-\bar{X})^2 = \frac{1}{n-1}(\sum_i X_i^2 -n\bar{X}^2)$）を用いると、

```math
\begin{equation}
    \mathcal{Z} \sim \int d\mu \left( \frac{n(\mu -\bar{X})^2}{(n-1)s^2} +1 \right)^{-\frac{n}{2}}
\end{equation}
```

となり、求めるt分布の表式を得ます。定数はこの積分が1になるように与えることができます。ここで通常$n=\nu+1$という変数を用います。（$\nu$は自由度と呼ばれます）

#F分布の導出

最後に、F分布の表式を導出します。F分布は2つの正規分布とその観測値$\{X_i\}$、$\{Y_i\}$が与えられたときの、標準偏差の比の分布です。
2つの正規分布$\{X_i\}$、$\{Y_i\}$の平均値をそれぞれ$\mu_x,\mu_y$、標準偏差をそれぞれ$\sigma_x,\sigma_y$と表します。さらに、観測値のサンプル数をそれぞれ$n_x,n_y$と表します。
ふたたびベイズの定理より、

```math
\begin{equation}
    P(\mu_x,\mu_y,\sigma_x,\sigma_y|{X_i},{Y_i}) \sim
    P({X_i},{Y_i}|\mu_x,\mu_y,\sigma_x,\sigma_y)P(\mu_x)P(\mu_y)P(\sigma_x)P(\sigma_y)
\end{equation}
```

となります。ここでも再び$\mu$と$\sigma$の無事前情報分布として$P(\mu) \sim const.$、$P(\sigma)\sim\frac{1}{\sigma}$を用いて、さらに$v=\frac{\sigma_x}{\sigma_y}$の値を固定した上でそれ以外の変数をすべて積分していきます。積分は、

```math
\begin{equation}
    \begin{split}
        \mathcal{Z} \sim & \int
        \mathop{ \underline{ d\mu_x \prod_i\frac{1}{\sqrt{2\pi\sigma_x^2}} 
        \exp\left(-\frac{(X_i-\mu_x)^2}{2\sigma_x^2} \right) } }
        \limits_{\text{①}}
        \mathop{ \underline{ d\mu_y \prod_i\frac{1}{\sqrt{2\pi\sigma_y^2}} 
        \exp\left(-\frac{(Y_i-\mu_y)^2}{2\sigma_y^2} \right) } }
        \limits_{\text{②}}
        \frac{d\sigma_x}{\sigma_x}\frac{d\sigma_y}{\sigma_y}
    \end{split}
\end{equation}
```

のように与えられます。
ここでまずは$\mu_x$,$\mu_y$から積分していきます。

```math
\begin{equation}
    \begin{split}
        ①=&
        \int \left( \frac{1}{\sqrt{2\pi\sigma_x^2}} \right)^{n_x} \exp\left( -\frac{n_x}{2\sigma_x^2} \left[ \left(\mu_x-\frac{1}{n_x}\sum_i X_i \right)^2 -\left(\frac{1}{n_x}\sum_i X_i \right)^2 + \frac{1}{n_x} \sum_i X_i^2 \right] \right) d\mu_x\\
        =& \left( \frac{1}{\sqrt{2\pi\sigma_x^2}} \right)^{n_x}\sqrt{\frac{2\pi\sigma_x^2}{n_x}}\exp\left( -\frac{n_x-1}{2\sigma_x^2} \left( \frac{1}{n_x-1}\sum_i X_i^2 -\frac{1}{n_x(n_x-1)}\left( \sum_i X_i\right)^2 \right)\right)\\
        =& \left( \frac{1}{\sqrt{2\pi\sigma_x^2}} \right)^{n_x-1}\frac{1}{\sqrt{n_x}} \exp\left( -\frac{n_x-1}{2\sigma_x^2}s_x^2 \right)
    \end{split}
\end{equation}
```

最後の行に至る式変形では、標本分散$s_x^2$の表式を用いました。同様に、

```math
\begin{equation}
\text{②}=
\left( \frac{1}{\sqrt{2\pi\sigma_y^2}} \right)^{n_y-1}\frac{1}{\sqrt{n_y}} \exp\left( -\frac{n_y-1}{2\sigma_y^2}s_y^2 \right)
\end{equation}
```
が得られます。
次に、$\sigma_x, \sigma_y$について、これらの比$v=\frac{\sigma_y}{\sigma_x}$を固定した上で積分を実行します。単純に$(\sigma_x, \sigma_y) \rightarrow (v, \sigma_y)$と変数変換をして、面積要素もそれに従って変換していくこともできますが、ここでは$\int \delta (v-\frac{\sigma_y}{\sigma_x}) dv$を元の表式に作用させ、$\sigma_x$について積分していくことにより計算していきます。（$\int \delta(f(x))dx = 1/f'(x)$を用います。）

```math
\begin{equation}
\begin{split}
    \mathcal{Z} \sim \int &
    \left( \frac{1}{\sqrt{2\pi\sigma_x^2}} \right)^{n_x-1}
    \left( \frac{1}{\sqrt{2\pi\sigma_y^2}} \right)^{n_y-1}
    n_x^{-\frac{1}{2}}n_y^{-\frac{1}{2}}\\
    &\exp\left( -\frac{n_x-1}{2\sigma_x^2} s_x^2 \right)
    \exp\left( -\frac{n_y-1}{2\sigma_y^2} s_y^2 \right)
    \left( -\frac{\sigma_x^2}{\sigma_y} \right)
    \frac{dv}{\sigma_x}
    \frac{d\sigma_y}{\sigma_y}
\end{split}
\end{equation}
```

さらに$\frac{1}{2\sigma_y^2}=\sigma_y'$と変数変換を行い、積分を実行します。

```math
\begin{equation}
    \begin{split}
        \mathcal{Z} \sim & \int (2\pi)^{-\frac{n_x+n_y}{2}-1}\sigma_y^{-n_x-n_y+2}
        v^{n_x-1} n_x^{-\frac{1}{2}} n_y^{-\frac{1}{2}}
        \exp\left( -\frac{n_x-1}{2\sigma_x^2}s_x^2 \right)
        \exp\left( -\frac{n_y-1}{2\sigma_y^2}s_y^2 \right)
        \left( -\frac{dv}{v}\frac{d\sigma_y}{\sigma_y} \right)\\
        \sim & \int (2\pi)^{-\frac{n_x+n_y}{2}-1}
        (2\sigma_y')^{\frac{n_x+n_y}{2}-2}v^{n_x-2}
         n_x^{-\frac{1}{2}} n_y^{-\frac{1}{2}}
         dv d\sigma_y'
         \exp\left( -v^2\sigma_y'(n_x-1) s_x^2 -\sigma_y'(n_y-1)s_y^2 \right)
    \end{split}
\end{equation}
```

さらに、$\sigma_y'$について積分を実行し、定数の係数を無視すると、

```math
\begin{equation}
\begin{split}
    \mathcal{Z} \sim & \int v^{n_x-2} \sigma_y'^{\frac{n_x+n_y}{2}-2}
    \exp\left( -\left[ v^2(n_x-1)s_x^2 +(n_y-1)s_y^2 \right]\sigma_y' \right) dvd\sigma'_y\\
    \sim & \int \Gamma\left( \frac{n_x+n_y}{2}-1 \right)
    \left[v^2(n_x-1)s_x^2 + (n_y-1)s_y^2 \right]^{-\frac{n_x+n_y}{2}+1}
    v^{n_x-2}dv
\end{split}
\end{equation}
```

が得られます。
さらに$w=v^2\frac{(n_x-1)s_x^2}{(n_y-1)s_y^2}$という表式を用いて式を簡単化すると、

```math
\begin{equation}
    \begin{split}
        \mathcal{Z} \sim & \int \left( w+ \frac{(n_y-1)s_y^2}{(n_x-1)s_x^2} \right)^{-\frac{n_x+n_y}{2}+1}w^{\frac{n_x}{2}-\frac{3}{2}}dw\\
        \sim & \int \left( w+ \frac{d_y s_y^2}{d_x s_x^2} \right)^{-\frac{d_x+d_y}{2}}w^{\frac{d_x}{2}-1}dw
    \end{split}
\end{equation}
```

となります。最終行に至る式変形では$d_x\equiv n_x-1$,$d_y\equiv n_y-1$という通常F分布を表すときに用いられるパラメーターをここでも用いました。これは観測値$X_i,Y_i$のそれぞれの自由度を表します。
よって、我々は約束した通りF分布の表式を求める事ができました。定数係数は、$\frac{d_y s_y^2}{d_xs_x^2} \equiv \alpha$という表式を用いた上で積分を実行することにより得られ、

```math
\begin{equation}
    \begin{split}
        \mathcal{Z} \sim & \int_0^\infty (x+\alpha)^{-d_x-d_y}w^{d_x-1}dw\\
        = & \left[ \frac{1}{(-d_x-d_y+1)}(w+\alpha)^{-d_x-d_y+1} w^{d_x-1} \right] - \int \frac{d_x-1}{(-d_x-d_y+1)}(w+\alpha)^{-d_x-d_y+1} w^{d_x-2}\\
        = & \cdots = \frac{(d_x-1)!}{(d_x+d_y-1)\cdots d_y}\alpha^{-d_y}
        = \frac{(d_x-1)!(d_y-1)!}{(d_x+d_y-1)!}\alpha^{-d_y}
    \end{split}
\end{equation}
```

となり、実際にF分布の係数も再現する事ができました。

#まとめ
最初に約束したように、我々は$\chi^2$分布、t分布、そしてF分布の表式をベイズの定理と無事前情報分布を用いて、不要な変数を積分することにより得ることができました。
母分布が正規分布であるとしたとき、$\chi^2$分布は母分布の標準偏差$\sigma$の分布となり、t分布は母分布の平均$\mu$の分布となり、F分布は2つの正規分布の標準偏差の比である$\sigma_x / \sigma_y$の分布となりました。
私の意見では、これらの分布や関連する検定がこのような見通しの良い方法で導かれること自体が重要であり、これらの分布に関する特性等を覚える必要は全くないだろうと思います。
このノートが、これらの分布の理解に費やす労力を減らすことを願います。

#参考文献
- https://en.wikipedia.org/wiki/Student%27s_t-distribution
- Gelman AB, Carlin JS, Rubin DB, et al. (1997). Bayesian Data Analysis (2nd ed.). Boca Raton: Chapman \& Hall. p. 68. ISBN 9780412039911.
- 松原望 統計学入門
- https://ja.wikipedia.org/wiki/%E3%82%B8%E3%82%A7%E3%83%95%E3%83%AA%E3%83%BC%E3%82%BA%E4%BA%8B%E5%89%8D%E5%88%86%E5%B8%83
- 統計Web 母分散の信頼区間の求め方1(χ2乗分布)　https://bellcurve.jp/statistics/course/9212.html
- 統計Web 標本とt分布 https://bellcurve.jp/statistics/course/8968.html
- 統計Web F分布表 https://bellcurve.jp/statistics/course/9932.html