Up

簡単な回帰分析

Visual C++プログラム例

 

変量yをq個の変量、・・・、の１次式

（１）                 

で表す回帰モデルについて考える。ここで、eは誤差項である。とおけば

                           

となり、を定数項（y-切片）を表すものとする回帰式となる。すなわち、式（１）は特殊な場合として定数項をもつ場合も表すことができる。

式（１）のモデルによる分析は、最小２乗法によるもの、最尤法によるもの、ベイズ的方法によるものがあげられるが、ここではガウス・マルコフのモデル（Gauss-Markoff setup）という枠組みで最小２乗法によって分析（Timm, 1975、Rao,1973）するプログラムを用意した。最尤法による分析の場合、データ数が十分に多いという条件の下で漸近定理が用いられる。しかしGauss-Markoff setupの最小２乗法の場合は、誤差項が独立に同じ正規分布に従うという仮定の下で、式（１）を行列で表したときの行列の階数の間の条件が満たされているときに正確に成り立つ定理に基づく分析であり、漸近定理によるものではない。行列の階数の間の条件は、ここではデータ数が独立変量の個数qより多いことを意味する。

式（１）のモデルにおいて、係数の値を求め、その帰無仮説「」の検定とp値、95％信頼区間を求めるプログラムSimpleRegression.exeを用意した。プログラムSimpleRegression.exeは実行形式のファイルなので、ダウンロード・実行することができる。ただし、.NET Framework 3.5対応である。プログラムのソースコードファイルなどをまとめた圧縮ファイルSimpleRegression.zipも用意したので、圧縮ファイル名SimpleRegression.zipをクリックしてダウンロードしたものを解凍すればVisual C++2008などで開くことができる。解凍は、ダウンロードしたファイルをマウス右ボタンでクリックして表示されるメニュから「展開」あるいは「解凍」などを選べば行うことができる。

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

図１

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

図２

 

　プログラムの入力データファイルは以下の形式で用意する。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

/

   2

                   I       X

    -0.045         1    -0.614

     4.142         1     1.033

    -0.699         1    -1.002

     2.386         1     0.490

     1.840         1    -0.016

     0.404         1    -0.782

     1.834         1     0.102

     2.358         1     0.383

     1.159         1    -0.104

     4.488         1     1.069

     1.956         1     0.330

     0.132         1    -0.889

     5.608         1     1.717

    -1.248         1    -0.552

     2.114         1     0.088

     1.071         1    -0.460

     1.940         1     0.087

     0.843         1    -0.132

     0.516         1    -0.630

     1.599         1     0.212

/

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

上のデータファイルは、

                           

の単回帰モデル用のものである。スラッシュ/で始まる行の次の行に独立変数の数qを書く。その次の行に、独立変数からまでのラベルを書く。上の例では、定数項を表すにＩ、独立変数xを表すにXがラベルとして書かれている。ラベルを並べた行の次の行からデータを

              y                      

の形式で並べて書く。上の例ではは定数項を表しているので常に１である。最後のデータの次の行は行の先頭をスラッシュで初めて、データの終了を示す。

　図２において、入力データファイル名を設定して「開く」ボタンをクリックすると、次に出力ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。

図３

 

適当なファイル名を設定する。既存のファイル名を設定すると、そのファイル名のファイルに上書きされる。ファイル名の設定後、「保存」ボタンをクリックすると計算が始まる。計算終了後、計算結果などがファイルに書き出されてファイルが保存された旨を表示するダイアログウィンドウが表示される（図４）。

図４

 

「OK」ボタンのクリックでプログラムの実行終了となる。プログラムの実行終了後、出力ファイルを開くと以下のようになっている（上の入力ファイル例の場合）。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

Data File = F:\yasuharu\www\LaCoocan\VCpp\simpleregression\SimpleRegression\SimpleRegression\DataX1.txt

 

Number of independent variables = 2

Number of cases = 20

 

Data...

                   I         X

    -0.045         1    -0.614

     4.142         1     1.033

    -0.699         1    -1.002

     2.386         1      0.49

      1.84         1    -0.016

     0.404         1    -0.782

     1.834         1     0.102

     2.358         1     0.383

     1.159         1    -0.104

     4.488         1     1.069

     1.956         1      0.33

     0.132         1    -0.889

     5.608         1     1.717

    -1.248         1    -0.552

     2.114         1     0.088

     1.071         1     -0.46

      1.94         1     0.087

     0.843         1    -0.132

     0.516         1     -0.63

     1.599         1     0.212

 

Qe = 4.83507

 

b(I) = 1.58203:          F(1,18) = 186.24243         p = 0.0000000   95%ConfIntvl = [1.33848, 1.82558]

b(X) = 2.29515:          F(1,18) = 185.19070         p = 0.0000000   95%ConfIntvl = [1.94082, 2.64948]

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

先頭に入力ファイル名が書き出されている。独立変数の個数とデータ数が書き出され、続いて読み込んだデータが読み込まれた変量のラベルの下に書き出される。続いて残差の平方和

                           

が書き出されている。この計算は、行列を扱うクラス型myMatを用いて次の簡単なコードで算出されている。

 

                 double Qe = (matTransp(y - X * b) * (y - X * b))[0][0];

 

　クラス型myMatは、ヘッダーファイルmymat.hにおいて宣言されている。その一部を本ウェブサイトの後の方にあげておく。ヘッダーファイルmymat.hは、圧縮ファイルSimpleRegression.zipに含まれている。

　出力ファイルの最後に、各係数の値がb(のラベル)の値として書き出され、続いて帰無仮説「」に対するＦ値、ｐ値、および95％信頼区間（岡本安晴「データ分析のための統計学入門」おうふう、2009、6.2節「区間推定」参照）が書き出されている。これらの計算は、最尤法における漸近定理に基づくものではなく、Gauss-Markoff setupの下でのexact probabilityである。

 

　本プログラムでは、各係数の値の帰無仮説「」に対する検定は、１つ１つ個別に行われている。そのため、多重共線性（multicollinearity）に関わる問題が生じることがある。次のモデル

                           

で生成したデータの場合、との相関が高いとき、との係数がいずれも帰無仮説「」を棄却できないということになりうる。圧縮ファイルにあるデータファイルData3a.txtはその例である。このファイルを入力ファイルとしてプログラムを実行すると次の結果を得る。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

b(y-int) = 53.47551:     F(1,46) = 428.01024         p = 0.0000000   95%ConfIntvl = [48.27257, 58.67845]

b(X1) = 1.28968:         F(1,46) = 2.28487           p = 0.1374812   95%ConfIntvl = [-0.42772, 3.00708]

b(X2) = 0.73587:         F(1,46) = 0.74944           p = 0.3911478   95%ConfIntvl = [-0.97515, 2.44688]

b(X3) = 0.87558:         F(1,46) = 102.55610         p = 0.0000000   95%ConfIntvl = [0.70155, 1.04962]

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

　いま、

                           

という変数x1x2を用意して（データファイル名Data3b.txt）分析すると以下の結果となる。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

b(y-int) = 53.47186:     F(1,47) = 436.27660         p = 0.0000000   95%ConfIntvl = [48.32175, 58.62197]

b(X1X2) = 2.02450:       F(1,47) = 777.19218         p = 0.0000000   95%ConfIntvl = [1.87841, 2.17060]

b(X3) = 0.87511:         F(1,47) = 104.46618         p = 0.0000000   95%ConfIntvl = [0.70287, 1.04736]

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

　変量x1x2の係数は有意に０でないことが示されている。

　上の例の場合、単純に平均値を用いたが、一般的には主成分分析などで相関の高い変量をまとめることが考えられる。

 

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝  mymat.h  ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

class myMat {

        public:

                int m, n;

                myVec * vv;

 

                myMat(): m(0), n(0), vv(NULL){  }

                myMat(int mm, int nn): vv(NULL) {

                        if (!((mm > 0) && (nn > 0))){

                                m = 0; n = 0; vv = NULL;

                        } else {

                                set_size(mm, nn);

                        }

                }

 

                myMat(const myMat & tmplt){

                        m = tmplt.m;

                        n = tmplt.n;

 

                        if (m > 0){

                                vv = new myVec[m];

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

                                        vv[i].v = new double[n];

                                        vv[i].vctrn = n;

                                       

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

                                                vv[i].v[j] = tmplt.vv[i].v[j];

                                        }

                                }

                        }

                        else vv = NULL;

                }

       

                virtual ~myMat(){ if (vv != NULL) delete[] vv; }

               

                //     (Re-)size the matrix

                int set_size( int mm, int nn );

               

 

                int get_size_row(){ return m; }

 

                int get_size_col(){ return n; }

 

                //       *this[]

                myVec & operator[](int mm);

 

                //       *this = R

                myMat & operator = (const myMat & R);

 

                //      *this = *this + R

                myMat & operator += (const myMat & R);

 

                //      *this = *this - R

                myMat & operator -= (const myMat & R);

 

                //      *this = *this * R

                myMat & operator *= (const myMat & R);

 

                //      *this = *this / R

                myMat & operator /= (const myMat & R);

};

 

 

//    Calculate L + R

myMat operator+ (const myMat & L, const myMat & R);

 

//    Calculate L - R

myMat operator- (const myMat & L, const myMat & R);

 

//    Calculate L * R

myMat operator* (const myMat & L, const myMat & R);

 

//    Calculate L / R

myMat operator/ (const myMat & L, const myMat & R);

 

//      Calculate a'

myMat matTransp( myMat & a );

 

//              Return myMat which is converted from A

myMat ConvToMat( double ** A, int m, int n );

 

//      myVecの書き出し

void WriteVec( myVec Vec );

 

//      myVecの書き出し：書き出し幅指定

void WriteVec( int w, myVec Vec );

 

//      myMatの書き出し

void WriteMat( myMat Mat );

 

//      myMatの書き出し：書き出し幅指定

void WriteMat( int w, myMat Mat );

 

//      myMatの読み込み

void ReadMat( myMat & Mat );

 

//    c = a + b

int matAdd( myMat & a, myMat & b, myMat & c, int m, int n );

 

 

//    逆行列の計算：　掃き出し法

//     a : 逆行列を求める行列；   invMat : 逆行列を返す行列

int calcInvMat( myMat a, myMat & invMat, int n, double rzero = 1.0e-15 );

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

 

参考文献

Rao, C. R. (1973)  Linear Statistical Inference and its Applications, 2nd ed.  John Wiley & Sons.

Timm, N. H. (1975) Multivariate Analysis with Applications in Education and Psychology. Brooks/Cole Publishing Company, Inc.

 

岡本安晴　作

Up