Up

平均と標準偏差のベイズ推定

 

平均 、分散 の正規分布 の母集団から 個のデータ、・・・、 が独立に採られたとする。

である。パラメータとして、分散の平方根である標準偏差 を使うことにする。このとき、独立な 個のデータ、・・・、 の得られる確率は

と書くことができる。 を と の関数と見たものは尤度関数と呼ばれ、 と書く。すなわち、

また、事後確率は

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

ただし、

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

なお、２項分布を例としたベイズ分析の説明を、＜岡本安晴「データ分析のための統計学入門」、pp. 160-167＞で行っている。

 

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

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

図１

 

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

図２

 

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

図３

 

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

図４

 

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

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

図５

 

左に並んでいる「Close」ボタンをクリックすると、図５のウィンドウは消えて、次に  の事後分布を示すグラフが表示される（図６）。

図６

 

左に並んでいる「Close」ボタンのクリックで図６のウィンドウは消えて、図７のメッセージボックスが表示され、図４で設定した出力用ファイルが保存されたことが示される（図７）。

図７

 

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

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

図８

 

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

 

 

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

 

提案分布

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

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

 

                              GenMuSgm *= scale_sgm(r_acpt_mean, n_param);

                        GenSDSgm *= scale_sgm(r_acpt_sd, n_param);

 

GenMuSgmとGenSDSgmが と であり、r_acpt_meanとr_acpt_sdは直前の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型で計算した場合、分子あるいは分母の確率が０になる可能性が高くなる。これは、double型の表しうる０でない数値の範囲が絶対値で 以上となっていることによる。すなわち、例えば、確率が0.1のデータを400個ぐらい独立に得たとき、その確率は となり、double型の範囲を超えた小さい値なので０として扱われる。このことを示すため、次のプログラムを用意した。

 

#include "stdafx.h"

 

using namespace System;

 

int main(array<System::String ^> ^args)

{

        double v = 1.0, Lv = 0.0;

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

                v *= 0.1;

                Console::WriteLine((i + 1).ToString()->PadRight(7) + "v = " + v.ToString());   

        }

        Console::ReadLine();

    return 0;

}

 

上のプログラムの実行結果を図９に示す。

図９

 

 を過ぎた行から精度が低下し、 を過ぎたあたりで０になっていることがわかる。Double型の有効桁数は１５である。この確率の積が０になることを避けるため、対数をとって計算している。すなわち、

である。このとき、

となる。

対数をとって計算した値を戻すときは指数関数を使うが、このときもdouble型の扱える数値の範囲による制約がある。これを調べるプログラムとして次のものを用意した。

 

#include "stdafx.h"

 

using namespace System;

 

int main(array<System::String ^> ^args)

{

        double v = Math::Round(Math::Log(1.0e305));

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

                v += 3.0;

                Console::WriteLine("Exp(" + v.ToString() + ") = " + Math::Exp(v)

                        + "   Exp(" + (-v).ToString() + ") = " + Math::Exp(-v));

        }

        Console::ReadLine();

    return 0;

}

 

上のプログラムの実行結果を図１０に示す。

図１０

 

 の の値が負の方向に絶対値が大きくなると、 の値は０になっているが、逆に の値が正の方向に絶対値が大きくなると の値は に設定されている。確率を計算しているので、 の値が１以上になるときは、 の値がいくら大きくなっても関数の値は１になる関数を用意して用いることにする。関数Exp0はそのための関数である。

 

        double Exp0(double v){

                if (v > 0.0){

                        return 1.0;

                }

                else{

                        return Math::Exp(v);

                }

        }

 

ステップ２における採択確率の計算もステップ１の場合と同様である。対数をとっての計算は以下のようになる。

ただし、

のときは、

となるので、上の対数による計算は行わず、採択確率を０とおいている。採択の判定は、一様乱数との比較で行っているので、比較の基準値 a の値を０以下の－１とおいて、実質的に採択確率を０としている。

 

上の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(mean[t], GenMuSgm);

                        double log_a = (-2.0 * sumX * (y - mean[t]) + n * (sqr(y) - sqr(mean[t])))

                                                                / (-2.0 * sqr(sd[t]));

                        double a = Exp0(log_a);

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

                                mean[t + 1] = y;

                                acpt_mean++;

                        }

                        else{

                                mean[t + 1] = mean[t];

                        }

 

                        y = rn->normalMS(sd[t], GenSDSgm);

                        if (y <= 0){

                                a = -1.0;

                        }

                        else{

                                log_a = n * Math::Log(sd[t] / y)

                                                        + 0.5 * ((1.0 / sqr(sd[t])) - (1.0 / sqr(y)))

                                                        * (ssumX - 2.0 * mean[t + 1] * sumX + n * sqr(mean[t + 1]));

                                a = Exp0(log_a);

                        }

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

                                sd[t + 1] = y;

                                acpt_sd++;

                        }

                        else{

                                sd[t + 1] = sd[t];

                        }

                }

                dmf.Close();

        }

 

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

提案分布  、 の調整は、 (GenMuSgm) および （GenSDSgm）の調整として、以下のコードで自動的に行っている。

 

                while (true){

                        MCMC ^ a_mcmc = gcnew MCMC();

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

                        a_mcmc->Init(nJump, n, X, NSimu, InitMean, InitSD, GenMuSgm, GenSDSgm);

                        a_mcmc->do_mcmc();

                        double sum_mean = 0.0, sum_sd = 0.0;

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

                                sum_mean += a_mcmc->mean[i + 1];

                                sum_sd += a_mcmc->sd[i + 1];

                        }

                        InitMean = sum_mean / NSimu;

                        InitSD = sum_sd / NSimu;

                        fout->WriteLine();

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

                        double r_acpt_mean = a_mcmc->acpt_mean / double(NSimu);

                        double r_acpt_sd = a_mcmc->acpt_sd / double(NSimu);

                        fout->WriteLine("r_acpt_mean = " + r_acpt_mean.ToString("F5") + "   r_accpt_sd = " + r_acpt_sd.ToString("F5"));

                        delete a_mcmc;

                        nJump++;

                        n_param = 0;

                        GenMuSgm *= scale_sgm(r_acpt_mean, n_param);

                        GenSDSgm *= scale_sgm(r_acpt_sd, 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