カイ２乗分布の計算
自由度dのカイ２乗分布（
-distribution）
に従う確率変数Xに対して
の関係にあるあるいはpを一方が与えられたとき他方を求めるプログラムcalc_chi2.exeを作成した。実行形式のプログラムcalc_chi2.exeは、.NET Framework3.5がインストールされているパソコンであれば、ファイル名calc_chi2.exeをクリックすることにより、ダウンロード、実行できる。プログラムのソースファイルなどをまとめたものは、圧縮ファイルcalc_chi2.zipとして用意したので、ファイル名calc_chi2.zipをマウスの右ボタンでクリックしてダウンロードできる。ダウンロードしたファイルはマウスの右ボタンで表示されるメニュから「解凍」あるいは「展開」を選んで解凍すれば、Visual C++などで開くことができる。なお、上の式ではガンマ関数を求めることになるが、ガンマ関数を求めずに計算するプログラムも＜ガンマ関数の計算を行わない方法＞に用意した。
プログラムは、実行すると図１のウィンドウ・フォームが表示される。
図１ 実行開始時のウィンドウ・フォーム
自由度dfとの値をそれぞれの右側のテキストボックスに設定して「chi2 -> p」ボタンをクリックすると、「p =」の右側のテキスト・ボックスに対応するpの値が表示される(図２)。
図２ 計算結果の表示。
「p =」の右側に値を設定して「p ->
xhi2」ボタンをクリックすると対応するの値が「Chi-Square =」の右側のテキスト・ボックスに表示される（図３）。
図３
累積確率は適応的数値積分によって求めているが、詳しくは
岡本安晴「統計学を学ぶための数学入門[下]」、2009、培風館、8.6節「積分の数値計算――ガウス・ルジャンドルの積分公式と適応的方法――」
などを参照されたい。
カイ２乗分布の式
により計算するプログラムを作成した。実行形式プログラムcalc_chi2a.exeは、Visual Studio 2013のインストールされているPCであれば、ファイル名calc_chi2a.exeをクリックして、ダウンロード・実行することがきる。プログラムのソースコードファイルなどをまとめた圧縮ファイルcalc_chi2a.zipも用意した。圧縮ファイル名calc_chi2a.zipをクリックしてダウンロード・解凍（展開）すればVisual Studio 2013で開くことができる。
なお、カイ２乗分布の分布関数は自由度をパラメータとしてもつので、分布関数を関数オブジェクトとして表し、自由度をそのメンバーとして設定し、積分などは関数テンプレートとして作成した。ソースコードは以下のようである。
//
// Yasuharu
Okamoto, 2015.12
//
#ifndef
CHISQRDISTRI_CK
#define
CHISQRDISTRI_OK
#include
"adapIntgrtn.h"
using
namespace adapInt;
namespace
ChiSqrD {
inline
double sqr(double v){
return
v * v;
}
class
Chi0 {
public:
double
df;
Chi0(){}
Chi0(double
vdf){
df
= vdf;
}
double
operator()(double chi){
double
v = ((df / 2.0) - 1.0) * Math::Log(chi) - (chi / 2.0);
return
Math::Exp(v);
}
};
class
Chi0t{
public:
double
df;
Chi0
chi0;
Chi0t(){}
Chi0t(double
vdf){
df
= vdf;
chi0.df
= df;
}
double
operator()(double t){
double
v = chi0(Math::Tan(t)) / sqr(Math::Cos(t));
return
v;
}
};
class
CumChi2 {
public:
double
df;
Chi0
chi0;
double
coeff;
CumChi2(){}
CumChi2(double
vdf){
df
= vdf;
chi0.df
= df;
Chi0t
chi0t(df);
TAdaptiveGL<Chi0t>
gl(chi0t);
double
IntgrlChi0 = gl.integral(0.0, Math::PI / 2.0);
coeff
= 1.0 / IntgrlChi0;
}
void
setCumChi2(double vdf){
df
= vdf;
chi0.df
= df;
Chi0t
chi0t(df);
TAdaptiveGL<Chi0t>
gl(chi0t);
double
IntgrlChi0 = gl.integral(0.0, Math::PI / 2.0);
coeff
= 1.0 / IntgrlChi0;
}
double
operator()(double chi2){
if
(df < 1.0) throw gcnew Exception("df = " + df.ToString() + "
< 1.0");
TAdaptiveGL<Chi0>
gl(chi0);
double
v = gl.integral(0.0, chi2);
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;
}
System::Diagnostics::Debug::WriteLine("L_b
= " + L_b.ToString() + "
m = " + m.ToString() + " U_b = " + U_b.ToString());
}
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);
}
// Newton法による根の計算
template<class
TF, class TDF>
void
TNewton(double & x, // 根の初期値、戻り値は根 f(x)=c
double
c, // 関数値
TF
& f, // 関数
TDF
& df, // 導関数
double
acc){
int
NStep = 0;
int
ck = 0;
double
x0;
do
{
NStep++;
if
(NStep > 20) { // 20回以上繰り返しのとき
ck
= 1;
break;
};
x0
= x;
x
= x0 + (c - f(x0)) / df(x0);
if
(x <= 0.0) x = 0.5 * x0;
if
(x > x0 * 2.0) x = x0 * 2.0;
}
while ((Math::Abs(x - x0) >= (acc * (Math::Abs(x) + Math::Abs(x0))))
&&
((Math::Abs(x)
+ Math::Abs(x0)) >= 1.0e-14)
);
if
((ck == 0) && (x > 1.0e-16)) return;
// Bisection法で計算
double
L_b = 1.0;
do
{
L_b
*= 0.5;
}
while (f(L_b) >= c);
double
U_b = 1.0;
do
{
U_b
+= 1.0;
}
while (f(U_b) <= c);
TBisection<TF>(f,
1.0, c, L_b, U_b, x, acc);
return;
}
class
Chi2_density{
public:
double
df;
double
coeff;
Chi0
chi0;
Chi2_density(){}
Chi2_density(double
vdf){
df
= vdf;
chi0.df
= df;
Chi0t
chi0t(df);
TAdaptiveGL<Chi0t>
gl(chi0t);
double
IntgrlChi0 = gl.integral(0.0, Math::PI / 2.0);
coeff
= 1.0 / IntgrlChi0;
}
void
set_df(double vdf){
df
= vdf;
chi0.df
= df;
Chi0t
chi0t(df);
TAdaptiveGL<Chi0t>
gl(chi0t);
double
IntgrlChi0 = gl.integral(0.0, Math::PI / 2.0);
coeff
= 1.0 / IntgrlChi0;
}
double
operator()(double chi2){
double
v = chi0(chi2);
return
v * coeff;
}
};
class
Calc_chi2_for_p {
public:
double
df;
CumChi2
cum_chi2;
Chi2_density
chi2_density;
Calc_chi2_for_p(){}
Calc_chi2_for_p(double
vdf){
df
= vdf;
cum_chi2.setCumChi2(df);
chi2_density.set_df(df);
}
double
operator()(double p){
double
x2 = 1.0;
TNewton<CumChi2,
Chi2_density>(x2, p, cum_chi2, chi2_density, 1.0e-12);
return
x2;
}
};
}
#endif