ティ分布の計算
自由度dのティ分布（t-distribution）
に従う確率変数Xに対して
の関係にあるtあるいはpを一方が与えられたとき他方を求めるプログラムcalc_t_distri.exeを作成した。実行形式のプログラムcalc_t_distri.exeは、.NET Framework3.5がインストールされているパソコンであれば、ファイル名calc_t_distri.exeをクリックすることにより、ダウンロード、実行できる。プログラムのソースファイルなどをまとめたものは、圧縮ファイルcalc_t_distri.zipとして用意したので、ファイル名calc_t_distri.zipをマウスの右ボタンでクリックしてダウンロードできる。ダウンロードしたファイルはマウスの右ボタンで表示されるメニュから「解凍」あるいは「展開」を選んで解凍すれば、Visual C++などで開くことができる。なお、自由度d の値が大きくて計算できないときは、本ウェブサイトの後半に用意したプログラムを試してみると計算できる可能性がある。
プログラムは、実行すると図１のウィンドウ・フォームが表示される。
図１ 実行開始時のウィンドウ・フォーム
自由度df（3.5など、整数でなくてもよい。）とtの値をそれぞれの右側のテキストボックスに設定して「t -> p」ボタンをクリックすると、「p =」の右側のテキスト・ボックスに対応するpの値が表示される(図２)。
図２ 計算結果の表示。
「p =」の右側に値を設定して「p -> t」ボタンをクリックすると対応するtの値が「t =」の右側のテキスト・ボックスに表示される（図３）。
図３
累積確率は適応的数値積分によって求めているが、詳しくは
岡本安晴「統計学を学ぶための数学入門[下]」、2009、培風館、8.6節「積分の数値計算――ガウス・ルジャンドルの積分公式と適応的方法――」
などを参照されたい。
ガンマ関数
を求めずに計算
自由度dのティ分布（t-distribution）を次式
の形で計算する。プログラムは、実行形式のファイルcalc_t_distri_1.exeとして用意したものは、ファイル名calc_t_distri_1.exeをクリックして実行することができる。使い方は、上の図１〜３の場合と同じである。プログラムのソースコードなどは、圧縮ファイルcalc_t_distri_1.zipとしてまとめたので、圧縮ファイル名calc_t_distri_1.zipをクリックダウンロードしたものを解凍すれば、Visual Studio 2013で開くことができる。上の式による計算では、かなり大きいd の値に対しても計算が可能である。自由度d の値が大きくなると正規分布に近づくが、これは正規分布のプログラムと比較することにより確認できる。
また、ティ分布の式は、パラメータとして自由度d をもつが、これを関数オブジェクトとして扱っている。すなわち、以下のようになっている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
// This
program was written in Visual C++ 2013
// by Yasuharu okamoto,
2015
#ifndef
CALC_T_CK
#define
CALC_T_CK
#include
"stdafx.h"
#include
"CalcIntgrt.h"
using
namespace myLib;
namespace
calc_t {
extern
double sqr(double v);
class
Ct0 {
public:
double
df;
double
operator()(double x){
double
v = 1.0 / (1.0 + sqr(x) / df);
return
Math::Pow(v, (df + 1.0) / 2.0);
}
};
class
Ct0_tan{
public:
double
df;
Ct0
t0;
double
operator()(double u){
return
t0(Math::Tan(u)) / sqr(Math::Cos(u));
}
};
class
CCum_t {
public:
double
df, coeff;
Ct0
t0;
Ct0_tan
t0_tan;
CCum_t(){}
CCum_t(double
vdf) : df(vdf) {
t0.df
= df;
t0_tan.df
= df;
t0_tan.t0.df
= df;
TAdaptiveGL<Ct0_tan>
gl(t0_tan);
double
v = gl.integral(-Math::PI / 2.0, Math::PI / 2.0);
coeff
= 1.0 / v;
}
void
set_df(double vdf){
df
= vdf;
t0.df
= df;
t0_tan.df
= df;
t0_tan.t0.df
= df;
TAdaptiveGL<Ct0_tan>
gl(t0_tan);
double
v = gl.integral(-Math::PI / 2.0, Math::PI / 2.0);
coeff
= 1.0 / v;
}
double
operator()(double t){
if
(t == 0.0) {
return
0.5;
}
else
{
TAdaptiveGL<Ct0>
gl(t0);
if
(t > 0.0) {
double
v = gl.integral(0.0, t);
v
= 0.5 + v * coeff;
return
v;
}
else
{
double
v = gl.integral(t, 0.0);
v
= 0.5 - v * coeff;
return
v;
};
}
}
};
// Bisection法による根の計算
template
<class TF>
void
Bisection(
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);
}
// Newton法による根の計算
template
<class TF, class TDF>
void
Newton(
double
& x, // 初期値
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);
}
while ((Math::Abs(x - x0) >= (acc * (Math::Abs(x) + Math::Abs(x0))))
&&
((Math::Abs(x)
+ Math::Abs(x0)) >= 1.0e-14)
);
if
(ck == 0) return;
// Bisection法で計算
double
L_b = 0.0;
do
{
L_b
-= 1.0;
}
while (f(L_b) >= c);
double
U_b = 0.0;
do
{
U_b
+= 1.0;
}
while (f(U_b) <= c);
Bisection(f,
1.0, c, L_b, U_b, x, acc);
return;
}
class
Ct_density {
public:
double
df;
double
coeff;
double
operator()(double t){
double
v = Math::Pow(1.0 / (1.0 + sqr(t) / df), (df + 1.0) / 2.0);
return
v * coeff;
}
};
class
Ccalc_t_for_p {
public:
double
df, coeff;
CCum_t
Cum_t_distribution;
Ct_density
t_density;
Ccalc_t_for_p(double
vdf) : df(vdf){
Cum_t_distribution.set_df(df);
coeff
= Cum_t_distribution.coeff;
t_density.df
= df;
t_density.coeff
= coeff;
}
double
operator()(double p){
double
t = 0.0;
Newton<CCum_t,
Ct_density>(t, p, Cum_t_distribution, t_density, 1.0e-12);
return
t;
}
};
}
#endif
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝