Up

ティ分布の計算

 

 自由度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 の値が大きくて計算できないときは、本ウェブサイトの後半に用意したプログラムを試してみると計算できる可能性がある。

 プログラムは、実行すると図1のウィンドウ・フォームが表示される。

図1 実行開始時のウィンドウ・フォーム

 

自由度df3.5など、整数でなくてもよい。)とtの値をそれぞれの右側のテキストボックスに設定して「t -> p」ボタンをクリックすると、「p =」の右側のテキスト・ボックスに対応するpの値が表示される(図2)

図2 計算結果の表示。

 

p =」の右側に値を設定して「p -> t」ボタンをクリックすると対応するtの値が「t =」の右側のテキスト・ボックスに表示される(図3)。

図3

 

累積確率は適応的数値積分によって求めているが、詳しくは

岡本安晴「統計学を学ぶための数学入門[]」、2009、培風館8.6節「積分の数値計算――ガウス・ルジャンドルの積分公式と適応的方法――」

などを参照されたい。

 

ガンマ関数を求めずに計算

 自由度dのティ分布(t-distribution)を次式

の形で計算する。プログラムは、実行形式のファイルcalc_t_distri_1.exeとして用意したものは、ファイル名calc_t_distri_1.exeをクリックして実行することができる。使い方は、上の図1〜3の場合と同じである。プログラムのソースコードなどは、圧縮ファイル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

=========================================================

 

 

Up