雑記: t分布の話

キャラクターは架空のものです。お気付きの点がありましたらご指摘いただけますと幸いです。
参考文献

  1. 24-3. 2標本t検定とは | 統計学の時間 | 統計WEB ― 2標本t検定の話があります。
  2. Student's t-distribution - Wikipedia ― t分布を発表された方の顔写真があります。
  3. http://www012.upp.so-net.ne.jp/doi/biostat/CT39/ttest.pdf ― 「t分布はどのようにして生まれたか」というのと「分散をプールするとはどういうことか」というのがわかりやすかったです。この記事のt分布の導入の流れはこの参考文献の流れに近いです。
  4. 意外と難しいT分布についての証明 - Qiita ― t分布の導出について全面的に参考にしました。
  5. 不偏分散と自由度n-1のカイ二乗分布 | 高校数学の美しい物語 ― 上の記事にもリンクがありますが不偏分散の分布の導出について全面的に参考にしました。

まとめ

  • 対応のない2標本(分散は同一で未知)の平均に差があるか検定するには自由度 n_0 + n_1 - 2 のt分布をつかう。
    • まず N(\mu_0, \sigma_0^2) から独立に n_0 個のデータを生成したら  (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがう。
    • 未知の \sigma_0^2 を不偏標本分散で置き換えた  (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} は自由度 n_0 - 1 のt分布にしたがう。
      •  (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがうことと、それとは独立に (n_0 - 1)S_0^2/\sigma_0^2 が自由度 n_0 - 1 のカイ2乗分布にしたがうことから導かれる。
    • 対応のない2標本があるとき、 \bigl( \bar{X_1} - \bar{X_0} - (\mu_1 - \mu_0) \bigr) / \sqrt{(1/n_0 + 1/n_1)\sigma_0^2}N(0,1) にしたがう。他方、(n_0 - 1)S_0^2/\sigma_0^2 + (n_1 - 1)S_1^2/\sigma_0^2 は自由度  n_0 + n_1 - 2 のカイ2乗分布にしたがう(カイ2乗分布の再生性)。これらを組み合わせると、自由度 n_0 + n_1 - 2 のt分布にしたがう検定統計量が構成される。

f:id:cookie-box:20200101101603p:plain:w60

日常生活の中で、手元にある対応のない2標本の平均に差があるか気になって、「標本0は N(\mu_0, \sigma_0^2) から独立に生成されていて、標本1は N(\mu_1, \sigma_0^2) から独立に生成されているだろう(\sigma_0 は不明)」と考えて、以下の検定統計量 T を用いて2標本t検定をすることがありますよね。但し、標本0のサイズを n_0, 標本平均を \bar{X_0}, 不偏標本分散を S_0^2 とし、標本1のサイズを n_1, 標本平均を \bar{X_1}, 不偏標本分散を S_1^2 とします。

 \displaystyle  T = \frac{\bar{X_1} - \bar{X_0}}{\sqrt{S^2 \left( \frac{1}{n_0} + \frac{1}{n_1}\right)}}, \quad S^2 = \frac{(n_0 - 1) S_0^2 + (n_1 - 1) S_1^2}{n_0 + n_1 - 2}
この T の実現値 t が棄却域にあるかどうかを調べればいいわけです。帰無仮説に応じた棄却域は以下です。ここで t_{a\%}(n) は自由度 n のt分布の a パーセント点です。
帰無仮説対立仮説棄却域(有意水準5%の場合)
 \mu_0 = \mu_1 \mu_0 \neq \mu_1 t < t_{2.5\%}(n_0 + n_1 - 2) または t_{97.5\%}(n_0 + n_1 - 2) < t
 \mu_0 \leqq \mu_1 \mu_0 > \mu_1t_{95\%}(n_0 + n_1 - 2) < t
自由度 n のt分布とは以下のような確率密度関数をもつ分布ですね。
 \displaystyle f(t) = \frac{\Gamma \left( \frac{n+1}{2} \right)}{\sqrt{n \pi} \,\Gamma \left( \frac{n}{2} \right)} \left( 1 + \frac{t^2}{n}\right)^{-\frac{n + 1}{2}}
でもよく考えると、どうしてこれで検定できるんでしょう? そもそもt分布というのもよく知らない分布ですし。どこからきた分布なんです?

f:id:cookie-box:20200101101614p:plain:w60

どこからきたというか、そこからきたんだよ。まず標本0にだけ注目しよう。いま標本0の母平均 \mu_0 の信頼区間を知りたいとする。標本0が N(\mu_0, \sigma_0^2) から独立に生成されていることより、 (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがうから、もし母分散 \sigma_0^2 がわかっていれば、適当な信頼係数で「 (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0} がとりうる区間」を決め打って、それを母平均について解けばいい。信頼係数 0.95 なら、信頼区間は標準正規分布の97.5%点 z_{97.5\%} を用いて以下のようになるね。

 \displaystyle \bar{X_0}- z_{97.5\%} \sqrt{\frac{\sigma_0^2}{n_0}} \leqq \mu_0 \leqq \bar{X_0}+ z_{97.5\%} \sqrt{\frac{\sigma_0^2}{n_0}}
でも母分散 \sigma_0^2 がわからなかったらどうする?

f:id:cookie-box:20200101101603p:plain:w60

母分散 \sigma_0^2 がわからなかったら、ですか…。先ほどと同様に  (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがうことを利用しようとしても、\sigma_0^2 がわからないので区間が定まりませんね。…であれば、\sigma_0^2 を不偏標本分散 S_0^2 などで置き換えて  (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} としてはいけないでしょうか? いえ、しかし、不偏標本分散 S_0^2 は標本の出方によって真の母分散 \sigma_0^2 よりも大きく出たり小さく出たりしますから、やはりそのまま置き換えても N(0, 1) にはしたがわないでしょう。不偏標本分散 S_0^2 が大きく出たり小さく出たりという分布をも考慮して、 (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} の分布の形を考え直す必要があるでしょうね…。

f:id:cookie-box:20200101101614p:plain:w60

それでいいよ。まず S_0^2 が大きく出たり小さく出たりという分布を出してみよう。

f:id:cookie-box:20200101101603p:plain:w60

えっ、S_0^2 のしたがう分布ですか? うーん…とりあえず展開しましょう。標本0内の i 番目のデータを X_{0, i} とかくことにしましょう( i = 1, \cdots, n_0 )。

 \displaystyle  \begin{split} S_0^2 &= \frac{1}{n_0 -1}\sum_{i=1}^{n_0} (X_{0, i} - \bar{X_0})^2 \\ &= \frac{1}{n_0-1}\sum_{i=1}^{n_0} X_{0, i}^2 -\frac{2 \bar{X_0}}{n_0 -1}\sum_{i=1}^{n_0} X_{0, i} + \frac{n_0}{n_0 -1} \bar{X_0}^2 \\ &= \frac{1}{n_0 -1}\sum_{i=1}^{n_0} X_{0, i}^2 - \frac{n_0}{n_0 -1} \bar{X_0}^2 \end{split}
第1項(第1項の中に n_0 項分あるわけですが)の分布と第2項の分布は2乗する前は正規分布ですから、変数変換すればこれらの分布も求ま…あ、第1項と第2項は独立ではないですね、負けました。

f:id:cookie-box:20200101101614p:plain:w60

急に投了しないで! …まあ確かに第2項は厄介だ。以下にもトリッキーってかいてあるしね。

つまり、座標変換で第2項を消すことにする。まず X_{0, i} がしたがう分布が N(\mu_0, \sigma_0^2) ではなくて N(0,1) の場合を考える。このとき Y_0 = Q X_0 という変換をする( X_0Y_0 は縦に n_0 個並んだベクトルのイメージね)。あ、この Q は直交行列にしたいのね。変換前後で2乗和を保ちたいから。じゃあどんな Q を選べばいいかというと、変換後の第1成分が  Y_{0,1} = (1/\sqrt{n_0}) \sum_{i=1}^{n_0}  X_{0, i} となるようにしたい。つまり、Q の1行目の成分をすべて 1/\sqrt{n_0} にする。だってそうすれば、
 \displaystyle  S_0^2 = \frac{1}{n_0 -1}\sum_{i=1}^{n_0} X_{0, i}^2 - \frac{n_0}{n_0 -1} \bar{X_0}^2 = \frac{1}{n_0 -1}\sum_{i=1}^{n_0} Y_{0, i}^2 - \frac{1}{n_0 -1} \sum_{i=1}^{n_0} \bar{X_{0,i}}^2 = \frac{1}{n_0 -1}\sum_{i=2}^{n_0} Y_{0, i}^2
と、厄介だった第2項が消えるからね。実際にこんな直交行列 Q を手に入れるには、単位行列の1行目の成分をすべて 1/\sqrt{n_0} にしてからグラムシュミットの直交化法をすればいいだろう。

f:id:cookie-box:20200101101603p:plain:w60

なんと…第1項は正規分布にしたがう確率変数の2乗和だから、直交変換をしても構わない。適当な直交変換を選べば、厄介な第2項が消える…。そして変換後の Y_0 = Q X_0 のしたがう分布はというと、多変量標準正規分布にしたがう確率ベクトルを直交変換しても多変量標準正規分布にしたがいますね。であれば (n_0 -1)S_0^2 = \sum_{i=2}^{n_0} Y_{0, i}^2 のしたがう分布は正規分布の二乗和がしたがう分布であり、つまり、自由度 n_0 -1 のカイ2乗分布ですか(これも導出すべきですが尺の都合で割愛します)。

f:id:cookie-box:20200101101614p:plain:w60

いまは X_{0, i} がしたがう分布が N(0,1) の場合を考えたけど、一般に X_{0, i}N(\mu_0, \sigma_0^2) にしたがう場合は (n_0 -1)S_0^2 / \sigma_0^2 が自由度 n_0 -1 のカイ2乗分布にしたがうね。

f:id:cookie-box:20200101101603p:plain:w60

…いま  (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} がしたがう分布を知ろうとして、S_0^2 のしたがう分布がわかった、という状態ですよね。では次はどうすればいいんでしょう?

f:id:cookie-box:20200101101614p:plain:w60

じゃあ分子の  \bar{X_0} のしたがう分布は?

f:id:cookie-box:20200101101603p:plain:w60

\bar{X_0} 自体は、 (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがうのでしょう? まず母分散 \sigma_0^2 が既知の場合はといって副部長自身が言っていましたよ。しかし、\bar{X_0}S_0^2 がしたがう分布がわかったからといって、 (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} のしたがう分布はわかりませんよね。結局、\bar{X_0} がどんな値をとるときは S_0^2 がどんな値をとりうるという相関構造がわからなければ…\bar{X_0}S_0^2 が独立でもない限り…。

f:id:cookie-box:20200101101614p:plain:w60

\bar{X_0}S_0^2 は独立だよ?

f:id:cookie-box:20200101101603p:plain:w60

え? 直感的には標本平均と標本分散が独立には思えな…。

f:id:cookie-box:20200101101614p:plain:w60

S_0^2 のしたがう分布を導出したときに、Y_0 = Q X_0 という変換をしたけど、変換後の Y_0 も多変量標準正規分布にしたがうから各成分は独立だ。そして、変換後の第1成分は  Y_{0,1} = (1/\sqrt{n_0}) \sum_{i=1}^{n_0} X_{0, i} = \sqrt{n_0} \bar{X_0} で、変換後の第2~n成分の2乗和は \sum_{i=2}^{n_0} Y_{0, i}^2 = (n_0 -1)S_0^2 だった。

f:id:cookie-box:20200101101603p:plain:w60

あっ…\bar{X_0}S_0^2 は独立ということになりますね…。Y_0 = Q X_0 という変換が標本平均と標本分散を振り分けていたとは…。であれば、いま知りたい  (\bar{X_0}-\mu_0) / \sqrt{S_0^2/n_0} のしたがう分布は計算すればよいです。

 \begin{split} \displaystyle P \left( \frac{\bar{X_0}-\mu_0}{\sqrt{S_0^2/n_0}} < t \right) &= P \left( \frac{\bar{X_0}-\mu_0}{\sqrt{\sigma_0^2/n_0}} \cdot \sqrt{\frac{\sigma_0^2}{(n_0-1) S_0^2}}< \frac{t}{\sqrt{n_0 - 1}} \right) \\ &= \iint_{\frac{a}{\sqrt{b}} < \frac{t}{\sqrt{n_0 - 1}} } z(a) \chi (b ; n_0 -1) da \, db \\ &= \iint_{\frac{a}{\sqrt{b}} < \frac{t}{\sqrt{n_0 - 1}} } \frac{e^{- a^2/2}}{\sqrt{2 \pi}} \frac{b^{(n_0 - 1)/2 - 1}e^{-b/2}}{2^{(n_0 - 1)/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} da \, db \end{split}
この変数を (a, b) から  (b, c=a/\sqrt{b}) に変換します。統計検定1級の教科書の1章でこれっぽい変数変換をちょいちょいやってるので何とかなるでしょう。この変換のヤコビアンは以下です。
 \begin{cases} a = h_0 (b, c) = \sqrt{b} c & \\ b = h_1 (b, c) = b & \end{cases}
 \left| \begin{array}{cc} \partial_b h_0(b, c) & \partial_c h_0(b, c) \\ \partial_b h_1(b, c) & \partial_c h_1(b, c) \end{array} \right| = \left| \begin{array}{cc} \frac{1}{2} b^{-1/2} c & \sqrt{b} \\ 1 & 0 \end{array} \right| = - \sqrt{b}
なので、このヤコビアンの絶対値をかければ変数変換できますね。
 \begin{split} \displaystyle P \left( \frac{\bar{X_0}-\mu_0}{\sqrt{S_0^2/n_0}} < t \right) &= \int_{b=0}^{\infty} \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \frac{e^{- b c^2/2}}{\sqrt{2 \pi}} \frac{b^{(n_0 - 1)/2 - 1}e^{-b/2}}{2^{(n_0 - 1)/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} \sqrt{b} \, db \, dc \\ &= \int_{b=0}^{\infty} \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \frac{b^{n_0/2 - 1}e^{- b c^2/2-b/2}}{ \sqrt{2 \pi} \cdot 2^{(n_0 - 1)/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} \, db \, dc \end{split}
なんだこれ…どうしろと…この記事を参照するとガンマ関数をつくるために、エクスポネンシャルの肩に注目して (b, c) から  (s=b c^2/2 + b/2, c) と変数変換するのですか。ではまたヤコビアンの計算です。
 \begin{cases} b = h_0 (s, c) = s/(c^2/2 + 1/2) & \\ c = h_1 (s, c) = c & \end{cases}
 \left| \begin{array}{cc} \partial_s h_0(s, c) & \partial_c h_0(s, c) \\ \partial_s h_1(s, c) & \partial_c h_1(s, c) \end{array} \right| = \left| \begin{array}{cc}  1/(c^2/2 + 1/2) & \ast \\ 0 & 1 \end{array} \right| = \displaystyle \left( \frac{c^2}{2} + \frac{1}{2} \right)^{-1}
これを適用すると…なるほど s積分でガンマ関数がつくれますね。
 \begin{split} \displaystyle P \left( \frac{\bar{X_0}-\mu_0}{\sqrt{S_0^2/n_0}} < t \right) &= \int_{s=0}^{\infty} \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \frac{s^{n_0/2 - 1}e^{- s}}{ \sqrt{2 \pi} \cdot 2^{(n_0 - 1)/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} \frac{ \left( \frac{c^2}{2} + \frac{1}{2} \right)^{-1} }{\left( \frac{c^2}{2} + \frac{1}{2} \right)^{n_0/2 - 1}}\, ds \, dc \\ &= \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \frac{ \Gamma \bigl( n_0/2 \bigr) }{ \sqrt{\pi} \cdot 2^{n_0/2} \Gamma \bigl( (n_0 - 1)/2 \bigr)} \left( \frac{c^2}{2} + \frac{1}{2} \right)^{ - n_0/2} \, dc \\ &= \frac{\Gamma \bigl( n_0/2 \bigr)}{\sqrt{\pi} \, \Gamma \bigl( (n_0 - 1)/2 \bigr)} \int_{c < \frac{t}{\sqrt{n_0 - 1}} } \left( c^2 + 1 \right)^{ - n_0/2 } \, dc \\ &= \frac{\Gamma \bigl( n_0/2 \bigr)}{\sqrt{(n_0 - 1) \pi} \, \Gamma \bigl( (n_0 - 1)/2 \bigr)} \int_{c' < t} \left( \frac{{c'}^2}{n_0 - 1} + 1 \right)^{ - n_0/2} \, dc' \end{split}
最後 c' = \sqrt{n_0 - 1} \, c とおきました。そしてこれは自由度 n_0 - 1 のt分布ですね…。

f:id:cookie-box:20200101101614p:plain:w60

ここまでたどり着いたことになるね。

定理1.標本0( サイズ n_0, 標本平均 \bar{X_0}, 不偏標本分散 S_0^2 )があるとき、以下の T_0 は自由度 n_0 - 1 のt分布にしたがう。
 \displaystyle T_0 = \frac{\bar{X_0} - \mu_0}{\sqrt{S_0^2 / n_0}}

f:id:cookie-box:20200101101603p:plain:w60

t分布が何者かはわかりました。まさしくその T_0 がしたがう分布なのですね。しかし、本当に検定したいのは標本0と標本1の平均が同じかどうかです。

f:id:cookie-box:20200101101614p:plain:w60

だったら調べるべき \mu_1 - \mu_0 をつくろう。いま以下がわかっていた。

  •  (\bar{X_0}-\mu_0) / \sqrt{\sigma_0^2/n_0}N(0,1) にしたがう。
  •  (\bar{X_1}-\mu_1) / \sqrt{\sigma_0^2/n_1}N(0,1) にしたがう。
ここから \mu_1 - \mu_0 をつくるには以下のように変形して正規分布の再生性をつかえばいい。
  •  \bar{X_0}N(\mu_0,\sigma_0^2/n_0) にしたがう。
  •  \bar{X_1}N(\mu_1,\sigma_0^2/n_1) にしたがう。
  •  \bar{X_1} - \bar{X_0}N(\mu_1 - \mu_0, \sigma_0^2/n_0 + \sigma_0^2/n_1) にしたがう。
  •  \bigl( \bar{X_1} - \bar{X_0} - (\mu_1 - \mu_0) \bigr) / \sqrt{(1/n_0 + 1/n_1)\sigma_0^2}N(0, 1) にしたがう。
もし \sigma_0^2 が既知なら、この信頼区間はただちに求まる。
 \displaystyle \bar{X_1} - \bar{X_0} - z_{97.5\%} \sqrt{ \left( \frac{1}{n_0} + \frac{1}{n_1} \right) \sigma_0^2} \leqq \mu_1 - \mu_0 \leqq \bar{X_1} - \bar{X_0} + z_{97.5\%} \sqrt{ \left( \frac{1}{n_0} + \frac{1}{n_1} \right) \sigma_0^2}
でも実際には \sigma_0^2 はわからない。だから S_0^2, S_1^2 で代用するしかない。いま上の議論より (n_0 - 1)S_0^2/\sigma_0^2 は自由度 n_0 - 1 のカイ2乗分布にしたがい、同様に (n_1 - 1)S_1^2/\sigma_0^2 も自由度 n_1 - 1 のカイ2乗分布にしたがう。カイ2乗分布もまた再生性をもつから、 (n_0 - 1)S_0^2/\sigma_0^2 + (n_1 - 1)S_1^2/\sigma_0^2 は自由度  n_0 + n_1 - 2 のカイ2乗分布にしたがう。このシチュエーションでさっきの議論と同様にt分布にたどり着くようにしたいなら以下のようにすればいい。今回はt分布の自由度も  n_0 + n_1 - 2 になるね。
 \displaystyle P \left( \frac{\bar{X_1} - \bar{X_0}-(\mu_1 - \mu_0)}{\sqrt{ \left( \frac{1}{n_0} + \frac{1}{n_1} \right) \sigma_0^2}} \cdot \sqrt{\frac{\sigma_0^2}{(n_0-1) S_0^2 + (n_1-1) S_1^2}}< \frac{t}{\sqrt{n_0 + n_1 - 2}} \right)
これを少し整理すればこうだ。帰無仮説 \mu_0 = \mu_1 のもとで2標本t検定の検定統計量に一致するよね。
 \displaystyle P \left( \frac{\bar{X_1} - \bar{X_0}-(\mu_1 - \mu_0)}{\sqrt{ \left( \frac{1}{n_0} + \frac{1}{n_1} \right) \frac{(n_0-1) S_0^2 + (n_1-1) S_1^2}{n_0 + n_1 - 2} }}< t \right)

つづかない