モチベーション
皆さんは, 非整数自由度の統計関数を使いたいと思ったことはありますか?
例えば, 正規分布モデルによる平均の区間推定を行う場合には, Studentのt分布を用いることになります. t分布は自由度$\nu$と変数$t$を入力として持っていますが, この自由度$\nu$は多くの場合, サンプルサイズ$-1$で自然数となります. カイ二乗分布やF分布の自由度に関しても, 多くの場合は自然数で事足ります.
しかし, Welchのt検定を行うとき, 非整数次の自由度を用いることになります. 統計学の教科書では「Welchのt検定では自由度を自然数に丸めてから計算する」と書いてある場合がありますが, 自由度を丸めてしまうと検定の棄却確率の精度が低下してしまいます. そもそも, t分布・カイ二乗分布・F分布の確率分布は, 正の数自由度について定義することができ, 自然数に制限する必要はありません.
多くの統計ソフトでは非整数自由度の統計関数を用いることができると思いますが, Excelの場合は(分析ツールを含む)標準関数ではそれができません. Excelにおけるt分布の確率密度関数や累積分布関数, その逆関数は, T.DISTやT.INVにより計算できますが, 自由度は整数で与えることしかできません.
この記事では, 非整数自由度のt分布, 加えてカイ二乗分布・F分布の計算をExcel等で計算する方法を示したいと思います. その根本原理は, それらの分布をベータ分布あるいはガンマ分布を用いて書き換えることにあります. Excelではベータ分布やガンマ分布は非整数のパラメータを受け付けるので, これが可能になるのです.
t分布
まずは, t分布をベータ分布に書き換えましょう.
自由度$\nu$のt分布の確率密度関数$f_\mathrm{t}(t;\nu)$は, 規格化定数を除いて,
\newcommand{\rt}{\mathrm{t}}
\newcommand{\id}{\mathrm{d}}
f_\rt(t;\nu)\propto \left(1+\frac{t^2}{\nu}\right)^{-\frac{\nu+1}{2}}
と書かれます. 今, $t\ge 0$としてt分布の対称性を考慮すると, 自由度$\nu$のt分布の(下側)累積分布関数$F_\rt(t;\nu)$は,
F_\rt(t;\nu)=\frac{1}{2}+\int_{0}^{t}\id u\,f_\rt(u;\nu)
と書くことができます.
$\int_{0}^{t}\id u f_\rt(u;\nu)$についてもう少し進めましょう.
\int_{0}^{t}\id u\,f_\rt(u;\nu) \propto \int_{0}^{t}\id u\,
\left(1+\frac{u^2}{\nu}\right)^{-\frac{\nu+1}{2}}
であり, $x=\left(1+\frac{u^2}{\nu}\right)^{-1}$と置換すると,
\newcommand{\rB}{\mathrm{B}}
\begin{split}
\int_{0}^{t}\id u\,f_\rt(u;\nu) &\propto
\int_{\frac{\nu}{\nu+t^2}}^{1}\id x\,
x^{\frac{\nu}{2}-1}(1-x)^{-\frac{1}{2}} \\
&=\rB\left(\frac{\nu}{2},\frac{1}{2}\right)
-\rB\left(\frac{\nu}{2},\frac{1}{2},\frac{\nu}{\nu+t^2}\right)
\end{split}
となります. ただし, $\rB(a,b)$はベータ関数, $\rB(a,b,z)$は不完全ベータ関数です.
さらに, $\int_{0}^{+\infty}\id u f_\rt(u;\nu) =\frac{1}{2}$となるように比例定数を定めると,
\int_{0}^{t}\id u\,f_\rt(u;\nu) =\frac{1}{2}-\frac{1}{2}
I\left(\frac{\nu}{2},\frac{1}{2},\frac{\nu}{\nu+t^2}\right)
となります. ただし, $I(a,b,z)=\rB(a,b,z)/\rB(a,b)$は正則化不完全ベータ関数です.
以上の議論により, $t\ge 0$のとき,
F_\rt(t;\nu)=1-\frac{1}{2}
I\left(\frac{\nu}{2},\frac{1}{2},\frac{\nu}{\nu+t^2}\right)
と表されることがわかりました. $t<0$のときはt分布の対称性により求めることができ, 合わせると
F_\rt(t;\nu)=\begin{cases}
\frac{1}{2}
I\left(\frac{\nu}{2},\frac{1}{2},\frac{\nu}{\nu+t^2}\right) &(t\le 0)\\
1-\frac{1}{2}
I\left(\frac{\nu}{2},\frac{1}{2},\frac{\nu}{\nu+t^2}\right) &(t\ge 0)
\end{cases}
となります.
正則化不完全ベータ関数$I(a,b,z)$はベータ分布の累積分布関数ですので, Excel上ではBETA.DIST(z,a,b,TRUE)関数を用いて実現されます.
次に, t分布の累積分布関数の逆関数を計算しましょう. 今, 与えられた$q$に対して$F_\rt(t;\nu)=q$となるような$t$を$t=F_\rt^{-1}(q;\nu)$と表します. $q\ge \frac{1}{2}$であれば$t\ge 0$ですので, まずは$q\ge \frac{1}{2}$として進めます. このとき,
I\left(\frac{\nu}{2},\frac{1}{2},\frac{\nu}{\nu+t^2}\right) =2-2q
が成立しますので, これを解いて
t=F_\rt^{-1}(q;\nu)=\sqrt{\nu\left(\frac{1}{I^{-1}\left(\frac{\nu}{2},\frac{1}{2},2-2q\right)}-1\right)}
が得られます. ただし, $I^{-1}(a,b,y)$は$I(a,b,z)$の逆関数:
\begin{gather}
z=I^{-1}(a,b,y) \\
\Longleftrightarrow\quad y=I(a,b,z)
\end{gather}
を表します.
$t<0$の場合も合わせて,
F_\rt^{-1}(q;\nu)=\begin{cases}
-\sqrt{\nu\left(\frac{1}{I^{-1}\left(\frac{\nu}{2},\frac{1}{2},2q\right)}-1\right)} &(q\le \frac{1}{2})\\
\sqrt{\nu\left(\frac{1}{I^{-1}\left(\frac{\nu}{2},\frac{1}{2},2-2q\right)}-1\right)} &(q\ge \frac{1}{2})
\end{cases}
となります.
正則化不完全ベータ関数の逆関数$I^{-1}(a,b,y)$はベータ分布の累積分布関数の逆関数ですので, Excel上ではBETA.INV(y,a,b)により実現されます. ただし, $q=\frac{1}{2}$のときはBETA.INVの計算結果が#NUM!になってしまうので, 条件分岐を加える必要があります.
最後に, t分布の確率密度関数$f_\rt(t;\nu)$を考えましょう. 累積分布関数$F_\rt(t;\nu)$:
F_\rt(t;\nu)=\begin{cases}
\frac{1}{2}
I\left(\frac{\nu}{2},\frac{1}{2},\frac{\nu}{\nu+t^2}\right) &(t\le 0)\\
1-\frac{1}{2}
I\left(\frac{\nu}{2},\frac{1}{2},\frac{\nu}{\nu+t^2}\right) &(t\ge 0)
\end{cases}
の両辺を$t$で微分することにより,
\newcommand{\pa}{\partial}
f_\rt(t;\nu)=
\frac{\pa F_\rt(t;\nu)}{\pa t}=
\frac{\nu|t|}{(\nu+t^2)^2}\left.\frac{\pa I(\frac{\nu}{2},\frac{1}{2},z)}{\pa z}\right|_{z=\frac{\nu}{\nu+t^2}}
が得られます. $\frac{\pa I(a,b,z)}{\pa z}$はベータ分布の確率密度関数ですので, Excel上ではBETA.DIST(z,a,b,FALSE)関数を用いて実現されます. ただし, $t=0$のときはBETA.DISTの計算結果が#NUM!になってしまうので, 条件分岐を加える必要があります. $f_\rt(0)$を陽に書くと
f_\rt(0;\nu)=\frac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\nu\pi}\Gamma(\frac{\nu}{2})}
です. ただし, $\Gamma(a)$はガンマ関数です.
$\Gamma(\frac{\nu+1}{2})$や$\Gamma(\frac{\nu}{2})$はそのまま計算するとオーバーフローしやすいので, 対数変換した
\frac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu}{2})}=\exp\left[
\textstyle \log\Gamma(\frac{\nu+1}{2}) - \log\Gamma(\frac{\nu}{2})
\right]
で計算することをお勧めします($\log$は自然対数).
対数ガンマ関数$\log\Gamma(a)$は, ExcelではGAMMALNで計算できます.
カイ二乗分布
自由度$\nu$のカイ二乗分布は, 形状母数(shape parameter)が$\frac{\nu}{2}$, 尺度母数(scale parameter)が$2$のガンマ分布に一致します. したがって, カイ二乗分布はGAMMA.DIST, GAMMA.INVにより実装することができます.
確率密度関数, 累積分布関数, 累積分布関数の逆関数はそれぞれ, GAMMA.DIST(x,nu/2,2,FALSE), GAMMA.DIST(x,nu/2,2,TRUE), GAMMA.INV(q,nu/2,2)で計算できます.
F分布
まずは, F分布はベータ分布を用いて書き換えましょう. 自由度$d_1,d_2$のF分布の確率密度関数$f_\mathrm{F}(x;d_1,d_2)$は規格化定数を除いて,
\newcommand{\rF}{\mathrm{F}}
f_\rF(x;d_1,d_2)\propto x^{-1}
\left(\frac{d_1 x}{d_1 x + d_2}\right)^{\frac{d_1}{2}}
\left(\frac{d_2}{d_1 x + d_2}\right)^{\frac{d_2}{2}}
と書かれます. 自由度$d_1,d_2$のF分布の下側累積分布関数$F_\rF(x;d_1,d_2)$は,
\begin{split}
F_\rF(x;d_1,d_2)&=\int_{0}^{x}\id u\, f_\rF(u;d_1,d_2) \\
&\propto \int_{0}^{x}\id u\,
u^{-1}\left(\frac{d_1 u}{d_1 u + d_2}\right)^{\frac{d_1}{2}}
\left(\frac{d_2}{d_1 u + d_2}\right)^{\frac{d_2}{2}} \\
&\propto \int_{0}^{\frac{d_1 x}{d_1 x + d_2}}\id v\,
v^{\frac{d_1}{2}-1}
(1-v)^{\frac{d_2}{2}-1}
\end{split}
と書かれます. ただし, $v=\frac{d_1 u}{d_1 u + d_2}$と置換しました.
したがって, 規格化定数を考慮すると,
F_\rF(x;d_1,d_2)=I( \textstyle \frac{d_1}{2},\frac{d_2}{2},
\frac{d_1 x}{d_1 x + d_2})
と書けることがわかります. $I$は正則化不完全ベータ関数です.
繰り返しになりますが, 正則化不完全ベータ関数$I(a,b,z)$はベータ分布の累積分布関数ですので, Excel上ではBETA.DIST(z,a,b,TRUE), その逆関数はBETA.INV(y,a,b)により実現されます.
最後に, F分布の確率密度関数$f_\rF(x;d_1,d_2)$を考えましょう. 累積分布関数$F_\rF(x;d_1,d_2)$の両辺を$x$で微分することにより,
f_\rF(x;d_1,d_2)=
\frac{\pa F_\rF(x;d_1,d_2)}{\pa x}=
\frac{d_1d_2}{(d_1 x + d_2)^2}\left.\frac{\pa I(\frac{d_1}{2},\frac{d_2}{2},z)}{\pa z}\right|_{z=\frac{d_1 x}{d_1 x + d_2}}
が得られます. $\frac{\pa I(a,b,z)}{\pa z}$はベータ分布の確率密度関数ですので, Excel上ではBETA.DIST(z,a,b,FALSE)関数を用いて実現されます.
おわりに
この記事では, 非整数自由度のt分布, カイ二乗分布, F分布をExcelで計算するための解決策の原理(理論)を示しました. Excelで実用的に使うためには, VBAプログラミングにより関数化することをお勧めします. この記事が少しでも皆さんの統計解析の助けになりますように.