U(u,g,N-r)分布に従うΛのp値
シミュレーションプログラム
分布に従う統計量
のp値を求めるプログラムSimuLambda.exeを作成した。プログラム名SimuLambda.exeをクリックしてダウンロード実行すると図１のフォームが表示される。
図１
３つのパラメータ値、u、g、N - r、と統計量の値を設定してGOボタンをクリックするとシミュレーションが始まる（図２）。
図２
シミュレーションの方法については、拙著「心理学データ分析と測定」の「３.Ａ 付録： Bartlett approximationの分布」を参照されたい。
分布に従うサンプルを生成するシミュレーションにおける試行数が「isimu = 」の値として表示されている。この値が30000試行に達するまでシミュレーションが続く。30000試行に達すると図３に示すフォームが表示される。
図３
中央に提示されているフォームは、図４に示すもので、計算中であることを示している。
図４
計算が終了すると図５のフォームとなる。
図５
分布は０の方向が棄却閾なので、設定したの値（図２、図５では0.5）より左側（小さい方）の確率
、シミュレーションではサンプルの比率、が
値となる。図５では、
であることが示されている。値の区間推定を２項分布によって行った95％信頼区間が、95%CIの値として示されている。なお、図５の曲線は、累積比率曲線であるので、の値に対するｐ値を示すグラフと見ることもできる。
図５の画面は、左上のSaveボタンをクリックするとファイルに保存することができる。
プログラムを圧縮ファイルSimuLambda.zipとしてまとめたものは、圧縮ファイル名SimuLambda.zipをクリックしてダウンロードすることができる。ダウンロードしたファイルはマウスの右ボタンでクリックすることにより解凍して、Visual C++ 2013で開くことができる。解凍したファイルを開いたとき、図６のようにエラーメッセージが表示されることがある。このときは、フォームのタブ「MyForm.h[デザイン]」の右端のX印をクリックして表示を消し、その後、ソリューションエクスプローラ内の「MyForm.h」のところをダブルクリックするとフォームが表示される。フォームが表示された後、プロジェクトをビルド実行することができる。
図６
プログラムのシミュレーションの部分を以下に示す。
private: System::Void buttonGO_Click(System::Object^ sender, System::EventArgs^ e) {
buttonGO->Enabled = false;
StatLambda = double::Parse(textBoxLambda->Text);
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);
if (calcInvMat(XX, InvXX, q) != 0){
MessageBox::Show("Error in calculating the inverse matrix os XX)");
return;
};
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);
if (calcInvMat(CXXC, InvCXXC, g) != 0){
MessageBox::Show("Error in calculating the inverse matrix of CXXC");
return;
};
myMat XB(N,p);
XB = X * B;
int NSimu = 30000;
ndata = NSimu;
Lambda = new double[NSimu];
chi2 = new double[NSimu];
int ck_simu;
for (int isimu = 0; isimu < NSimu; isimu++){
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;
double numV = DetMat(Qea), denV = DetMat(Qta);
if (numV <= 0.0) {
MessageBox::Show("N-r = " + Nr.ToString() + " の値が小さ過ぎます。");
buttonGO->Enabled = true;
return;
}
Lambda[isimu] = numV / denV;
}
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();
}