エフ分布の計算
自由度(df1,df2)のエフ分布（F distribution）
に従う確率変数Xに対して
の関係にあるFとpにおいて、一方が与えられたとき他方を求めるプログラムcalc_F_distri.exeを作成した。実行形式のプログラムcalc_F_Distri.exeは、.NET Framework3.5がインストールされているパソコンであれば、ファイル名calc_F_distri.exeをクリックすることにより、ダウンロード、実行できる。プログラムのソースファイルなどをまとめたものは、圧縮ファイルcalc_F_distri.zipとして用意したので、ファイル名calc_F_distri.zipをマウスの右ボタンでクリックしてダウンロードできる。ダウンロードしたファイルはマウスの右ボタンによるクリックで表示されるメニュから「解凍」あるいは「展開」を選んで解凍すれば、Visual C++ 2008などで開くことができる。なお、このプログラムで扱える自由度より大きい値の自由度に対応できるプログラムをガンマ関数を直接計算しない方法で作成し、この後の＜ガンマ関数の計算を行わない方法＞に用意した。
プログラムは、実行すると図１のウィンドウ・フォームが表示される。
図１ 実行開始時のウィンドウ・フォーム
自由度df1とdf2、およびFの値をそれぞれの右側のテキストボックスに設定して「F -> p」ボタンをクリックすると、「p =」の右側のテキスト・ボックスに対応するpの値が表示される(図２)。自由度は、整数値でなくてもよいことに注意。
図２ 計算結果の表示。
「p =」の右側に値を設定して「p -> F」ボタンをクリックすると対応する
の値が「F =」の右側のテキスト・ボックスに表示される（図３）。
図３
累積確率は適応的数値積分によって求めているが、詳しくは
岡本安晴「統計学を学ぶための数学入門[下]」、2009、培風館、8.6節「積分の数値計算――ガウス・ルジャンドルの積分公式と適応的方法――」
などを参照されたい。
F分布を
によって計算するプログラムを作成した。実行形式のプログラムcalc_F_distri1.exeは、Visual Studio 2013がインストールされているPCであれば、プログラム名calc_F_distri1.exeをクリックしてダウンロード・実行できる。ソースプログラムは圧縮ファイルcalc_F_distri1.zipとして用意した。圧縮ファイル名calc_F_distri1.zipをクリックしてダウンロード解凍（展開）すれば、Visual Studio 2013で開くことができる。本プログラムでは自由度df1が１の場合は対象としていない。自由度df1が１の場合は、自由度df2のｔ分布の２乗が自由度（１、 df2）のＦ分布になることを利用してｔ分布のプログラムを利用すれば確率を求めることができる。
プログラムでは、関数オブジェクトとテンプレートを用いて、F分布のパラメータ（自由度）をオブジェクト内に取り込んで扱った。ソースコードは以下のようである。
//
// Y. Okamoto, 2015.12
//
#ifndef CALCFDISTRI_CK
#define CALCFDISTRI_CK
#include "adapIntgrtn.h"
using namespace adapInt;
namespace calcfd {
inline double sqr(double v){
return v * v;
}
class CF0 {
public:
double df1, df2;
CF0(){}
CF0(double vdf1, double vdf2){
df1 = vdf1;
df2 = vdf2;
if (df1 < 2) throw gcnew Exception("df1 should be at least 2.");
}
void SetDF(double vdf1, double vdf2){
df1 = vdf1;
df2 = vdf2;
if (df1 < 2) throw gcnew Exception("df1 should be at least 2.");
}
double operator()(double x){
if (df1 < 2) throw gcnew Exception("df1 should be at least 2.");
double v = ((df1 / 2.0) - 1.0) * Math::Log(x)
- ((df1 + df2) / 2.0) * Math::Log(1.0 + (df1 * x / df2));
return Math::Exp(v);
}
};
class CF0tan {
public:
double df1, df2;
CF0 f0;
CF0tan(){}
CF0tan(double vdf1, double vdf2){
df1 = vdf1;
df2 = vdf2;
f0.df1 = df1;
f0.df2 = df2;
}
double operator()(double t){
double v = f0(Math::Tan(t)) / sqr(Math::Cos(t));
return v;
}
};
class CCum_F {
public:
double df1, df2;
CF0 f0;
double coeff;
CCum_F(){}
CCum_F(double vdf1, double vdf2){
df1 = vdf1;
df2 = vdf2;
f0.df1 = df1;
f0.df2 = df2;
CF0tan f0tan(df1, df2);
TAdaptiveGL<CF0tan> gl(f0tan);
double vg0 = gl.integral(0.0, Math::PI / 2.0);
coeff = 1.0 / vg0;
}
double operator()(double F){
TAdaptiveGL<CF0> gl(f0);
double v = gl.integral(0.0, F);
return v * coeff;
}
};
// Bisection法による根の計算
template <class TF>
void TBisection(
TF & f, // 関数オブジェクト
double s, // fが増加関数のとき s > 0.0
double c, // 根に対する関数値
double L_b, double U_b, // 根の存在区間の下限と上限
double & Root, // 根 f(Root) = c
double acc){
if ((f(L_b) - c) * (f(U_b) - c) > 0.0){
throw gcnew Exception("上限と下限の組み合わせが不適切です。");
}
do {
double m = 0.5 * (L_b + U_b);
double v = f(m);
if (s * (v - c) > 0.0) {
U_b = m;
}
else {
L_b = m;
}
} while ((Math::Abs(U_b - L_b) >= acc * Math::Abs(U_b))
&&
((Math::Abs(U_b) + Math::Abs(L_b)) >= 1.0e-14)
);
Root = 0.5 * (L_b + U_b);
}
}
#endif