Up

２変量の平均値、標準偏差と相関係数

 

平均、分散、相関係数 の２変量正規分布の母集団から個のデータ、・・・、 が独立に採られたとする。すなわち、

である。、、、、および をパラメータとして確率密度関数を表せば

となる。したがって、個の独立なデータ、・・・、 が得られる確率は

となる。

 をパラメータ、、、、および の関数とみたものは尤度関数と呼ばれ、 と書く。

である。

このとき、事後確率は

で与えられる。データ数が多くなると、事後確率は主に尤度関数によって決まるようになる。すなわち、データの情報は尤度関数によって表される。データ数が多くないときも、事後確率が尤度関数によって決まるようにするため、事前確率を定数と置く。すなわち、

ただし、

とおき、 を一様分布とする。ただし、 が０でない領域（サポート）にわたる積分が１になるようにするため、サポートは十分に広い有界な領域であるとする。コンピュータでの乱数は有界な範囲での値しかとらないので、その範囲をカバーする十分な大きさのサポートであるとすれば、実質的には何ら制約を受けないことになる。本ウェブサイトのプログラムで用いている乱数についてはこのウェブサイトで説明している。

上の事後確率を求めるプログラムを作成した。プログラムは圧縮ファイルNormalTwoVars.zipとしてまとめたので、ファイル名NormalTwoVars.zipをクリックしてダウンロード解凍すれば、Visual Studio 2013で開くことができる（Visual Studio 2015でも可；Visual StudioはMicrosoft社のウェブサイトから無料でダウンロード インストールできる）。Visual Studio 2013のインストールされているPCであれば、実行形式のプログラムファイルNormalTwoVars.exeをダウンロード実行することもできる。Windows 8.1において実行形式のファイルをダウンロード実行するときの手順についてはこのウェブサイトで説明している。

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

図１

 

「GO」ボタンのクリックで、入力データファイル名の設定を求めるダイアログ ウィンドウが表示される（図２）。

図２

 

入力データファイルは、図３のようにスラッシュ / で始まる行の次の行から、１行に１組ずつデータを書いて用意する。最後のデータの次の行の先頭文字はスラッシュ / としてデータの終わりを示す。

図３

 

入力データファイル名の設定後、「開く」ボタンをクリックすると、次に計算結果などが出力される出力用ファイルの名前の設定を求めるダイアログ ウィンドウが表示される（図４）。

図４

 

適当なテキストファイル名を設定して、「保存」ボタンをクリックすると、計算が始まる。

計算が終了すると、まず の事後分布を示すグラフが表示される（図５）。

図５

 

左に並んでいる「Close」ボタンをクリックすると、図６のフォームが表示される。

図６

 

ボタン「rho」、「mu1」、「sigma1」、「mu2」、「sigma2」をクリックすると、それぞれ、、、、 の事後分布のグラフが表示される。例えば、「rho」ボタンをクリックすると図７の画面になる。

図７

 

図６のフォームにおいて、「Close」ボタンをクリックすると、図８のメッセージボックスが表示される。

図８

 

「OK」ボタンのクリックでプログラムは終了する。

プログラムの終了後、出力ファイルを開くと図９のようになっている。

図９

 

事後分布の平均値、４分位数、モードなどが出力されている。

 

プログラムでは、式（１）で与えられている事後確率を、Metropolis-within-GibbsによるMCMC（Markov chain Monte Carlo）によって求めている。Metropolis-within-Gibbsについては、岡本安晴「心理学データ分析と測定」勁草書房の付録C「ベイズ的分析法」pp. 227-232で説明している。式（１）の場合は、以下の手順になる。

 

まず、提案分布

を用意する。提案分布は、ウェブサイト＜ベイズデータ分析＞で説明されている形では、例えば は、 と なども条件に含まれる形であるが、は などに関しては定数であるものを選ぶので条件には表示していない。

　 などの値は、MCMCにおける更新確率（acceptance probability）が適切な範囲になるように調節する。具体的には、関数scale_sgmにより、以下のように適応的に調整している。

 

                                          GenMuSgm1 *= scale_sgm(r_acpt_mean1, n_param);

                        GenMuSgm2 *= scale_sgm(r_acpt_mean2, n_param);

                        GenSDSgm1 *= scale_sgm(r_acpt_sd1, n_param);

                        GenSDSgm2 *= scale_sgm(r_acpt_sd2, n_param);

                        GenRhoSgm *= scale_sgm(r_acpt_rho, n_param);

 

GenMuSgm1などが などであり、r_acpt_mean1などは直前のMCMCにおける更新比率である。関数scale_sgmは以下のように与えられているので、更新比率が0.3～0.5の範囲になるように自動的に調整される。

 

                                           double scale_sgm(double acpt_r, int & ck){

                                 double v = 1.0;

 

                                 if (acpt_r > 0.5){

                                         v = 1.0 + 9.0 * (acpt_r - 0.5) / 0.5;

                                         ck++;

                                 }

                                 else if (acpt_r < 0.3){

                                         v = 1.0 / (1.0 + 9.0 * (0.3 - acpt_r) / 0.3);

                                         ck++;

                                 }

 

                                 return v;

                         }

 

上の提案分布を用いて、以下のステップを繰り返す。

 

ステップ０

               とおき、初期値、、、 および を設定する。最初のMCMCにおける初期値は、データの平均値および標準偏差を初期値として設定する。提案分布における などの値も、データの標準偏差を基に設定する。２回目以降のMCMCにおいては、直前のMCMCでの平均値および調整値を用いる。

ステップ１

               の更新候補値を

ととる。

採択確率 を次式

により設定する。

              確率 で、

              確率 で

とおく。

ステップ２

               の更新候補値を

ととる。

採択確率 を次式

により設定する。

              確率 で、

              確率 で

とおく。

ステップ３

               の更新候補値を

ととる。

採択確率 を次式

により設定する。

              確率 で、

              確率 で

とおく。

ステップ４

               の更新候補値を

ととる。

採択確率 を次式

により設定する。

              確率 で、

              確率 で

とおく。

ステップ５

               の更新候補値を

ととる。

採択確率 を次式

により設定する。

              確率 で、

              確率 で

とおく。

ステップ６

             

とおく。 の値が目標数に達していたら終了するが、達していないときはステップ１に戻る。

 

ここで、上のアルゴリズムにおける計算上の注意をいくつか上げておく。詳しくは、ウェブサイト＜平均と標準偏差のベイズ推定＞で説明している。

ステップ１における採択確率の計算において

の部分は、

であることから、

である。ステップ２～５においても同様である。

また、

の部分は、データ が多くなるとC++のdouble型で計算した場合、分子あるいは分母の確率が０になる可能性が高くなる。この確率の積が０になることを避けるため、対数をとって計算している。すなわち、

である。このとき、

ここで、

、、、、

である。

ステップ３～５の場合、パラメータの値の範囲に制約がおかれているが、更新候補値が範囲外の値のときは、選択確率が０になるようにしている。

 

上のMCMCのアルゴリズムを表すコードは以下のようになる。

 

        void MCMC::do_mcmc(){

                dmf.Show();

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

                        if ((t % 1000) == 0){

                                dmf.ShowMessage(" " + t.ToString() + "/" + NSimu.ToString(), " ");

                        }

 

                        double y = rn->normalMS(mean1[t], GenMuSgm1);

                        double num = LogL(y, mean2[t], sd1[t], sd2[t], rho[t]);

                        double den = LogL(mean1[t], mean2[t], sd1[t], sd2[t], rho[t]);

                        double a = Exp0(num - den);

                        if (rn->uni() < a){

                                mean1[t + 1] = y;

                                acpt_mean1++;

                        }

                        else {

                                mean1[t + 1] = mean1[t];

                        }

 

                        y = rn->normalMS(mean2[t], GenMuSgm2);

                        num = LogL(mean1[t + 1], y, sd1[t], sd2[t], rho[t]);

                        den = LogL(mean1[t + 1], mean2[t], sd1[t], sd2[t], rho[t]);

                        a = Exp0(num - den);

                        if (rn->uni() < a){

                                mean2[t + 1] = y;

                                acpt_mean2++;

                        }

                        else {

                                mean2[t + 1] = mean2[t];

                        }

 

                        y = rn->normalMS(sd1[t], GenSDSgm1);

                        if (y <= 0.0){

                                a = 0.0;

                        }

                        else {

                                num = LogL(mean1[t + 1], mean2[t + 1], y, sd2[t], rho[t]);

                                den = LogL(mean1[t + 1], mean2[t + 1], sd1[t], sd2[t], rho[t]);

                                a = Exp0(num - den);

                        }

                        if (rn->uni() < a){

                                sd1[t + 1] = y;

                                acpt_sd1++;

                        }

                        else {

                                sd1[t + 1] = sd1[t];

                        }

 

                        y = rn->normalMS(sd2[t], GenSDSgm2);

                        if (y <= 0.0) {

                                a = 0.0;

                        }

                        else {

                                num = LogL(mean1[t + 1], mean2[t + 1], sd1[t + 1], y, rho[t]);

                                den = LogL(mean1[t + 1], mean2[t + 1], sd1[t + 1], sd2[t], rho[t]);

                                a = Exp0(num - den);

                        }

                        if (rn->uni() < a){

                                sd2[t + 1] = y;

                                acpt_sd2++;

                        }

                        else{

                                sd2[t + 1] = sd2[t];

                        }

 

                        y = rn->normalMS(rho[t], GenRhoSgm);

                        if ((y <= -1.0) || (y >= 1.0)) {

                                a = 0.0;

                        }

                        else{

                                num = LogL(mean1[t + 1], mean2[t + 1], sd1[t + 1], sd2[t + 1], y);

                                den = LogL(mean1[t + 1], mean2[t + 1], sd1[t + 1], sd2[t + 1], rho[t]);

                                a = Exp0(num - den);

                        }

                        if (rn->uni() < a){

                                rho[t + 1] = y;

                                acpt_rho++;

                        }

                        else{

                                rho[t + 1] = rho[t];

                        }

                }

 

上の関数do_mcmc は、参照クラス型MCMCのメソッドとして与えられている。参照クラス型については、＜岡本安晴「大学生のための心理学VC++プログラミング入門」勁草書房、第１１章 参照型＞で説明している。

提案分布 などの調整は、（GenMuSgm1）などの調整として以下のコードで自動的に行っている。

 

                while (true){

                        MCMC ^ a_mcmc = gcnew MCMC();

                        a_mcmc->dmf.SetMsg("Adaptive stage-" + (++n_step).ToString());

                        a_mcmc->Init(nJump, n, X, NSimu, InitMean1, InitMean2, InitSD1, InitSD2, InitRho,

                                GenMuSgm1, GenMuSgm2, GenSDSgm1, GenSDSgm2, GenRhoSgm);

                        a_mcmc->do_mcmc();

 

                        double sum_mean1, sum_mean2, sum_sd1, sum_sd2, sum_rho;

                        sum_mean1 = sum_mean2 = sum_sd1 = sum_sd2 = sum_rho = 0.0;

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

                                sum_mean1 += a_mcmc->mean1[i + 1];

                                sum_mean2 += a_mcmc->mean2[i + 1];

                                sum_sd1 += a_mcmc->sd1[i + 1];

                                sum_sd2 += a_mcmc->sd2[i + 1];

                                sum_rho += a_mcmc->rho[i + 1];

                        }

                        InitMean1 = sum_mean1 / NSimu;

                        InitMean2 = sum_mean2 / NSimu;

                        InitSD1 = sum_sd1 / NSimu;

                        InitSD2 = sum_sd2 / NSimu;

                        InitRho = sum_rho / NSimu;

                        fout->WriteLine();

                        fout->WriteLine("mean1 = " + InitMean1.ToString("F5") + "   sd1 = " + InitSD1.ToString("F5"));

                        fout->WriteLine("mean2 = " + InitMean2.ToString("F5") + "   sd2 = " + InitSD2.ToString("F5"));

                        fout->WriteLine("rho = " + InitRho.ToString("F5"));

 

                        double r_acpt_mean1 = a_mcmc->acpt_mean1 / double(NSimu);

                        double r_acpt_mean2 = a_mcmc->acpt_mean2 / double(NSimu);

                        double r_acpt_sd1 = a_mcmc->acpt_sd1 / double(NSimu);

                        double r_acpt_sd2 = a_mcmc->acpt_sd2 / double(NSimu);

                        double r_acpt_rho = a_mcmc->acpt_rho / double(NSimu);

                        fout->WriteLine();

                        fout->WriteLine("r_acpt_mean1 = " + r_acpt_mean1.ToString("F5"));

                        fout->WriteLine("r_acpt_mean2 = " + r_acpt_mean2.ToString("F5"));

                        fout->WriteLine("r_acpt_sd1 = " + r_acpt_sd1.ToString("F5"));

                        fout->WriteLine("r_acpt_sd2 = " + r_acpt_sd2.ToString("F5"));

                        fout->WriteLine("r_acpt_rho = " + r_acpt_rho.ToString("F5"));

 

                        delete a_mcmc;

                       

                        nJump++;

                        n_param = 0;

                        GenMuSgm1 *= scale_sgm(r_acpt_mean1, n_param);

                        GenMuSgm2 *= scale_sgm(r_acpt_mean2, n_param);

                        GenSDSgm1 *= scale_sgm(r_acpt_sd1, n_param);

                        GenSDSgm2 *= scale_sgm(r_acpt_sd2, n_param);

                        GenRhoSgm *= scale_sgm(r_acpt_rho, n_param);

                        fout->WriteLine("n_param = " + n_param.ToString());

                        if (n_param == 0) break;

                }

 

提案分布と初期値は、適応的に自動化された手続きで設定されるので、気にせずにプログラムを実行して分析を行うことができる。

MCMCによるサンプリングの結果は、カーネル スムーザーによって滑らかな曲線が求められる。これは次の関数k_smootherによって行われている。

 

        double k_smoother(double x, double * samples, int N, double sgm){

                double v = 0.0;

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

                        if (Math::Abs(x - samples[i]) < 4 * sgm){

                                v += Math::Exp(-0.5 * sqr((x - samples[i]) / sgm));

                        }

                return v / (N * sgm * Math::Sqrt(2.0 * Math::Atan(1.0) * 4.0));

        }

 

スムージングのときの幅は、パラメータdouble sgm によって与えられている。この値が小さいとグラフはサンプルの頻度の変化に敏感になる。例えば、図５において、kwの値（この値がパラメータsgmの値として使われる）を0.1に設定して（図１０）、「ReDraw」ボタンをクリックすると図１１のようになる。

 

図１０

 

図１１

 

逆に大きい値、例えば１０を設定して、「ReDraw」ボタンをクリックすると図１２のようになる。

図１２

 

グラフに表されている曲線からモードが算出されているので、山の頂上付近が滑らかな曲線になる範囲で、できるだけ小さい値をkw（関数k_smootherのパラメータsgmの値として使われる）に設定するようにする。kw の値を小さくすると、グラフはサンプルの頻度の変化に敏感になり精度はよくなるが、激しく変化するグラフになる（図１１）。kw の値を大きくすると、グラフは滑らかになるが、広い範囲での平均値となるので、精度が悪くなる（図１２）。

設定されたkwの値ごとに算出されたモードの値は、出力ファイルに出力されている（図１３）。

図１３

 

 

Up