Biometrics Bulletin 1946 [2] を最初に読んだけど、詳しい話が書いてなかったので Psychometrika 1941 [1] を読んだ。
概要は、不偏分散の推定量
\[
V_i=\frac{1}{\nu_i}\sum_{j=1}^{n_i}(Y_{ij}-\bar{Y}_i)^2
\] の和の分布および、カイ二乗分布に従う確率変数の線形結合の分布を近似する方法である。
[1] に書かれた数式の展開がよく分からなかったため、自分で色々と考えた結果をメモします。そのうち修正するかも。
まず、最初に $V_i$ の従う分布を示している (以降しばらく $i=1, 2$ としておく)。
ここで、
\[
X_i=\nu_i\frac{V_i}{\s_i^2}=\sum_{j=1}^{n_i}\left(\frac{Y_{ij}-\bar{Y}_i}{\s_i}\right)^2\sim \chi^2(\nu_i)
\] \[
\mathrm{d}x_i=\frac{\nu_i}{\s_i^2}{\mathrm d}v_i
\]
とおくと、$V_i$ の従う分布は
\[
\begin{eqnarray}
f(V_i=v_i)&=&
\frac{2^{-\nu_i/2}}{\Gamma(\nu_i/2)}\left(\nu_i\frac{v_i}{\s_i^2}\right)^{\nu_i/2-1}\exp\left(-\nu_i\frac{v_i}{2\s_i^2}\right)\frac{\nu_i}{\s_i^2}\\ &=&
\frac{1}{\Gamma(\nu_i/2)}\left(\frac{\nu_i}{2\s_i^2}\right)^{\nu_i/2}v_i^{\nu_i/2-1}\exp\left(-\nu_i\frac{v_i}{2\s_i^2}\right)\end{eqnarray}
\] になる。
これより、$Z=V_1+V_2$ の従う分布を考えると、
\[
\begin{eqnarray}
f(Z=z)&=&
\int_{0}^{+\infty}f_1(V_1=v_1)f_2(Z-V_1=z-v_1){\mathrm d}v_1 \\ &=&
K \exp\left( -2b\frac{z}{2\s_2^2} \right)\int_{0}^{+\infty}v_1^{a-1}(z-v_1)^{b-1}\exp\left( -2a\frac{v_1}{2\s_1^2}-2b\frac{-v_1}{2\s_2^2} \right){\mathrm d}v_1
\end{eqnarray}
\] が得られる。
ただし、
\[
K=\frac{1}{\Gamma(\nu_1/2)\Gamma(\nu_2/2)}\left(\frac{\nu_1}{2\s_1^2}\right)^{\nu_1/2}\left(\frac{\nu_2}{2\s_2^2}\right)^{\nu_2/2}
\] は規格化定数?であり、$a=\nu_1/2 \geq b=\nu_2/2$ とする。
上式は、$\nu_2$ が偶数の場合、$(z-v_1)^{b-1}$ を二項定理により展開することができ
\[
\begin{eqnarray}
f(Z=z)&=&
K \exp\left( -2b\frac{z}{2\s_2^2} \right) \sum_{k=0}^{b-1} z^k\begin{pmatrix}b-1\\k\end{pmatrix} \int_{0}^{+\infty}v_1^{(a-1)-(b-1)+k}\exp\left\{ -2(a\s_2^2-b\s_1^2)\frac{v_1}{2\s_1^2\s_2^2} \right\}{\mathrm d}v_1
\end{eqnarray}
\] となる。
さらに、$\nu_1$が偶数の場合、
\[
\begin{eqnarray}
f(Z=z)&=&
K \exp\left( -2b\frac{z}{2\s_2^2} \right) \sum_{k=0}^{b-1} z^k\begin{pmatrix}b-1\\k\end{pmatrix}\frac{\Gamma(a-b+k+1)}{\left(\frac{a}{\s_1^2}+\frac{b}{\s_2^2}\right)^{a-b+k+1}} \\ &=&
\sum_{k=0}^{b-1} K \begin{pmatrix}b-1\\k\end{pmatrix}\frac{\Gamma(a-b+k+1)}{\left(\frac{a}{\s_1^2}+\frac{b}{\s_2^2}\right)^{a-b+k+1}} (\s_2^2/2b)^k (2bz/\s_2^2)^{(k+1)-1} \exp\left( -\frac{2bz/\s_2^2}{2} \right)
\end{eqnarray}
\] と書けるので、おそらく自由度 $k+1$ のカイ二乗分布を混合したような複雑な分布になる事が分かる。
実際使うとなると、$\nu_1, \nu_2$ が偶数だろうが奇数だろうがこの分布は複雑すぎる。
だから、自由度 $\nu$ のカイ二乗分布で近似する事にする (近似できそうだという事を数値例を出して、グラフを描いて確かめている[1])。
具体的には、近似する分布の分散と理論的な $Z$ の分散が一致するように、$\nu$ を定める方針にする。
ここでは、
\[
Z=k V=\sum_{i=1}^2 V_i
\] とした。
右辺は不偏分散の和になっているので、左辺も同様に不偏分散の様なもの $V$ と何らかの係数 $k$ で表現できるものと考える事にする。
ここで、$V$ について、$V_i$ において $\nu_i, \s_i^2$ に対応する値を、それぞれ $\nu, \s^2$ と置くことにする。
$X_i$ の期待値と分散が、
\[
\mathrm{E}(X_i)=\nu_i,\quad \mathrm{Var}(X_i)=2\nu_i
\] だから、$V_i$ の期待値と分散は、
\[
\mathrm{E}(V_i)=\mathrm{E}\left(\frac{X_i\s_i^2}{\nu_i}\right)=\s_i^2
\] \[
\mathrm{Var}(V_i)=\mathrm{Var}\left(\frac{X_i\s_i^2}{\nu_i}\right)=\frac{2\s_i^4}{\nu_i}
\] となる。
これの類似性から、$Z$ の期待値と分散は、
\[
\mathrm{E}(Z)=\mathrm{E}(k V)=k \s^2, \qquad
\mathrm{Var}(Z)=\mathrm{Var}(k V)=\frac{2 (k \s^2)^2}{\nu}
\] と書く事ができるし、また、定義より、
\[
\mathrm{E}(Z)=\mathrm{E}\left(\sum_{i=1}^2 V_i\right)=\sum_{i=1}^2 \s_i^2, \qquad
\mathrm{Var}(Z)=\mathrm{Var}\left(\sum_{i=1}^2 V_i\right)=2\sum_{i=1}^2\frac{\s_i^4}{\nu_i}
\] と書く事もできる。
この4つの式を連立すると、
\[
\nu=\frac{\left(\sum_{i=1}^2 \s_i^2\right)^2}{\sum_{i=1}^2\frac{\s_i^4}{\nu_i}}
\] が得られる。
これが Satterthwaite の近似である。
和ではなく、線形結合
\[
Z=\sum_{i=1}^{n} k_i V_i
\] の場合にも、同様の手順で、
\[
\nu=\frac{\left(\sum_{i=1}^n k_i \s_i^2\right)^2}{\sum_{i=1}^n k_i^2 \frac{\s_i^4}{\nu_i}}
\] を導くことができる。
実際用いるときは、$\s_i^2$ は未知なので、不偏推定量で置き換えるなどの対応をする。
多分、一般的な状況において、このような形をした nuisance parameter をうまく消す方法はない。
