Up

U分布の バートレット カイ2乗近似

 

分布をで近似するバートレット近似（Bartlett’s approximation）

の精度をシミュレーションで調べるプログラムSimuBartlettUChi2.exeを作成した。プログラム名SimuBartlettUChi2.exeをクリックして、ダウンロード実行することができる。Windows8/8.1で表示されるセキュリティ確認メッセージについてはウェブサイト＜http://y-okamoto-psy1949.la.coocan.jp/VCpp/Win8EXEFile/＞で説明している。

プログラムの実行開始時に「***.dllがインストールされていません」というようなエラーメッセージが表示されたときは、Visual Studio Express 2013 for Windows Desktopをマイクロソフトのウェブサイト

http://www.visualstudio.com/downloads/download-visual-studio-vs

からダウンロードインストールすれば必要なファイルがインストールされ、実行できるようになる。

プログラムを実行すると図１のフォームが表示される。

図１

 

u、g、N – r に値を設定して「GO」ボタンをクリックするとシミュレーションが始まる（図２）。N は、多変量重回帰モデルにおけるデータ数（ケース数）であり、r はデザイン行列Xの階数であるが、U分布においてはN – r の値がパラメータの１つである。シミュレーションの方法については、拙著「心理学データ分析と測定」の「３.Ａ　付録： Bartlett approximationの分布」を参照されたい。

図２

 

シミュレーションにおける試行数が「isimu = 」の値として表示される。30000は総試行数である。isimuの値が30000に達すると図３のフォームが表示される。

図３

 

シミュレーションによって生成されたバートレット近似の値の累積頻度曲線が破線で描かれ、対応する分布の累積確率分布関数が実線で描かれている。図３の場合は、バートレット近似の分布は分布でよく近似されていることがわかる。

図の上部に表示されている

Prop(samples >= Chi2C95) = 4.88%　および　Prop(samples >= Chi2C99) = 0.91%　は、分布の95%と99%の値

に対して、シミュレーションにおけるサンプル値がそれぞれの値以上であったものの比率である。図３の場合、それぞれ５％および１％に近い値である。

図３の画面は、左上の「Save」ボタンをクリックすることによって、ファイルに保存することができる。保存したファイルは、ペイントなどで図４のように表示することができる。

図４

 

図４のようにu とg の値が大きくなると、データ数N の値も大きくとる必要があることがわかる。既に集めたデータの分析の場合は、データ数をバートレット近似が使える数まで増やすことは現実的ではない。本ウェブサイトのプログラム（ソースコードファイルなどをまとめたものは、このウェブサイトの後半に用意してある）を少し変更すれば、シミュレーションで直接分布の p 値を求めることができる。分布の p 値を求めるプログラム例のウェブサイトも用意した。

また、N – r の値が小さいとき、シミュレーションの過程で多変量重回帰分析の適切な解が求まらないことがあるが、その場合は、N – r の値を大きくするようにメッセージが表示される。

なお、u と g の値も、積 u*g が200以下という制限を設けているが、これは、分布の累積確率を数値計算する上での制限である。累積確率を本プログラムで用いているものとは異なる方法で算出する場合は、この制限とは異なる値となる。シミュレーションで直接分布の p 値を求める場合は、この制限はない。

 

多変量重回帰分析は、以下の手順で解を求めている。

多変量重回帰モデルを

とおく。ここで、Y は N 行 p 列の行列であり、はN 行 q 列の行列である。また、

とおく。

このとき、解を

とおくと

である。このときの誤差を

とおく。

に対するコントラストを

とおき、

とおく。

コントラストの制約条件のもとでの誤差を与える解

に対して

とおく。このとき、次式

が成り立ち、

である。

　これらの式に基づいて、シミュレーションプログラムは以下のようになっている。

 

         u = int::Parse(textBoxU->Text);       

         g = int::Parse(textBoxG->Text);       

         int Nr = int::Parse(textBoxNr->Text);

 

         if ((u * g) > 200){

                 MessageBox::Show("u * g = " + (u*g).ToString() + " の値が200以下になるようにして下さい。");

                 buttonGO->Enabled = true;

                 return;

         }

 

         p = u + 1;

         q = g + 1;

         r = q;

 

         N = Nr + r;   

 

         int df = u * g;

 

         myMat B(q,p);

         for (int i = 0; i < q; i++)

                 for (int j = 0; j < p; j++)

                         B[i][j] = 1.0;

 

         RNP2to191 rn;

 

         myMat X(N,q);

         for (int i = 0; i < N; i++)

                 for (int j = 0; j < q; j++)

                         X[i][j] = rn.uni() * 10.0;

 

         myMat A(p, u);

         for (int i = 0; i < p; i++)

                 for (int j = 0; j < u; j++){

                         if(i == j){

                                 A[i][j] = 1.0;

                         }else if( (i-1) == j){

                                 A[i][j] = -1.0;

                         }else{

                                 A[i][j] = 0.0;

                         }

                 }

 

        myMat C(g, q);

        for (int i = 0; i < g; i++)

                for (int j = 0; j < q; j++){

                        if (i == j){

                                C[i][j] = 1.0;

                        }else if((i + 1) == j){

                                C[i][j] = -1.0;

                        }else{

                                C[i][j] = 0.0;

                        }

                }

 

        myMat XX(q,q);

        XX = matTransp(X) * X;

        myMat InvXX(q,q);

        calcInvMat( XX, InvXX, q );

        myMat XXX(q, N);

        XXX = InvXX * matTransp(X);

        myMat XXXX(N,N);

        XXXX = X * XXX;

        myMat I_XXXX(N,N);

        for (int i = 0; i < N; i++)

                for (int j = 0; j < N; j++)

                        if(i == j){

                                I_XXXX[i][j] = 1.0 - XXXX[i][j];

                        }else{

                                I_XXXX[i][j] = -XXXX[i][j];

                        }

 

        myMat CXXC(g,g);

        CXXC = C * InvXX * matTransp(C);

        myMat InvCXXC(g,g);

        calcInvMat( CXXC, InvCXXC, g );

 

        myMat XB(N,p);

        XB = X * B;

 

        int NSimu = 30000;

 

        ndata = NSimu;

 

        double * Lambda = new double[NSimu];

        chi2 = new double[NSimu];

 

        int ck_simu;

        for (int isimu = 0; isimu < NSimu; isimu++){

                ck_simu = 0;

                do{

                        if ((isimu % 100) == 0){

                                labelMsg->Text = "isimu = " + isimu.ToString() + "/" + NSimu.ToString();

                                Update();

                        }

                        else if ((N > 500) && ((isimu % 10) == 0)){

                                labelMsg->Text = "isimu = " + isimu.ToString() + "/" + NSimu.ToString();

                                Update();

                        }

 

                        myMat Y(N,p);

                        for (int i = 0; i < N; i++)

                                for (int j = 0; j < p; j++)

                                        Y[i][j] = XB[i][j] + rn.normal();

 

                        myMat Bhat(q, p);

                        Bhat = XXX * Y;

                                       

                         myMat YA(N,u);

                         YA = Y * A;

                         myMat Qea(u,u);

                         Qea = matTransp(YA) * I_XXXX * YA;

 

                         myMat CBhatA(g,u);

                         CBhatA = C * Bhat * A;

                         myMat Qha(u,u);

                         Qha = matTransp(CBhatA) * InvCXXC * CBhatA;

                         myMat Qta(u,u);

                         Qta = Qha + Qea;

 

                         Lambda[isimu] = DetMat(Qea) / DetMat(Qta);

 

                         chi2[isimu] = -((N - r) - 0.5 * (u - g + 1)) * Math::Log(Lambda[isimu]);

 

                         if((0 < chi2[isimu]) && (chi2[isimu] < 1.0e4)){

                                 ck_simu = 1;

                         }else{

                                 MessageBox::Show("N-r の値が小さ過ぎます。");

                                 buttonGO->Enabled = true;

                                 return;

                         };

                }while(ck_simu == 0);

        }

 

        double sum = 0.0;

        for (int i = 0; i < NSimu; i++)

                sum += chi2[i];

        double meanChi2 = sum / NSimu;

 

         DrawGraph ^ dg = gcnew DrawGraph();

         dg->SetValues( NSimu, Lambda, u * g );

         dg->pf = this;

         dg->Show();

         dg->Draw();

 

プログラムのソースコードファイルなどを圧縮ファイルSimuBartlettUChi2.zipにまとめた。ファイル名SimuBartlettUChi2.zipをクリックしてダウンロード解凍すれば、Visual Studio 2013で開くことができる。

ダウンロード解凍したファイルをVisual Studio 2013で開いたとき、図５のようにエラーが表示されることがある。

図５

 

このエラーメッセージは無視してよい。ソリューションエクスプローラの「MyForm.h」（図６）をダブルクリックすれば、フォームが表示される。F5キーを押して、ビルド実行することもできる。

図６

 

 

Up