多次元正規分布に従う乱数列を生成する方法

[質問]

多次元正規分布に従う乱数列を生成するにはどのようにしたらよいでしょうか。

[回答]

2次元までであれば、DATAステップだけで比較的簡単に作成することができます。

変数 X および Y の母平均をμx, μy、母標準偏差をσx, σy、分散共分散行列をσxyとします。このようなパラメータの正規分布は、次のような回帰モデルに基づいて、平均0分散1の独立な2つの確率変数ε1およびε2から生成することができます。

  Ｘ＝μx + σx * ε1
  Ｙ＝μy + a* (x-μx) + b * ε2

上式において、a=σxy/σx^2 , b=sqrt(σy^2 - σxy^2/ σx^2)です。このような式に基づいて、下記のようなプログラムで乱数を作成することができます。

  DATA mnormal0;
    KEEP x y;

    /* パラメータ */
    mux=2;  muy=1;          /* 母平均 */
    varx=4; vary=5;         /* 母分散 */
    covxy=2;                /* 母共分散 */
    n=100;                  /* 標本の大きさ */
    seed=1234;              /* シード値 */

    /* エラー処理 */
    IF (vary-covxy**2/varx)<=0 THEN DO;
      PUT 'ERROR:分散共分散行列が正定値(positive definite)ではありません。';
      STOP;
    END;

    /* 乱数生成のためのループ */
    DO i=1 TO n;
      x = mux+SQRT(varx)*RANNOR(seed);
      y = muy + covxy / varx* (x-mux) + SQRT(vary -covxy**2/varx)*RANNOR(seed);
      OUTPUT;
    END;
  RUN;

  /* 生成された乱数列の確認 */
  PROC CORR DATA=mnormal0 COV;
    VAR x y;
  RUN;

一般に、多変量正規分布に従う乱数列の生成は、分散共分散行列をCholesky分解することに基づいて行うことができます。
このCholesky分解は、SAS/STATのMIXEDプロシジャに備わる機能を応用して実現可能です。

次の例は、3次元正規分布に従う乱数列を100オブザベーション作成しています。なお、Mersenne-Twister (MT)による乱数生成は、SAS 8.2では評価版の機能です。

                                                      /* MTの例 */

  %LET nobs=100;                          /* オブザベーション数 */
  %LET ncol=3;                                          /* 次元 */
  %LET seed=12345;                                    /* シード */
  %LET out=mnormal1;                      /* 出力データセット名 */

             /* 平均と共分散行列からなる SAS データセットを作成 */
                       /* 1列目から3列目までは共分散行列の要素  */
                              /* 平均は最後の列で「縦」に入れる */
  DATA cov1;
    INPUT col1-col3 mean;
    Row+1;
  DATALINES;
   1    1.2 2.25 20
   1.2  4   3.3  30
   2.25 3.3 9    40
  ;
  run;

                 /* 標準正規分布に従う乱数を100×3 だけ作成する */
                  /* 関数 RAND ではなく関数 RANNOR を使用すると */
              /* 乗算型合同法に基づく乱数列を生成することになる */
  DATA _random;
    CALL STREAMINIT(&seed);
    _TYPE_="SCORE";
    _MODEL_="col";
    mean=1;
    ARRAY col{&ncol};
    DO num=1 TO &nobs;
      DO i=1 TO DIM(col);
        col{i}=RAND("NORMAL");
      END;
      output;
    END;
    DROP i;
  RUN;

                                       /* MIXEDプロシジャの利用 */
                    /* MIXEDプロシジャは、Cholesky 分解を行なう */
        /* 目的のためにのみ使用しており、プログラムそのものには */
                                /* 統計的な意味づけはできません */
  ODS LISTING CLOSE;
  ODS OUTPUT CHOLG=_Cholesky;
  PROC MIXED DATA=cov1;
    CLASS Row mean;
    PARMS /NOITER;
    MODEL mean=;
    RANDOM row*mean/TYPE=UN GDATA=cov1 GC; 
  RUN;
  ODS LISTING;


                /* SCORE プロシジャで行列の乗算を擬似的に行なう */
  PROC SCORE DATA=_cholesky SCORE=_random OUT=_out(KEEP=col num);
    BY num;
    VAR mean col:;
  RUN;

  PROC TRANSPOSE DATA=_out OUT=&out.(DROP=_NAME_);
    BY num;
  RUN;

なお、SAS/ETSのMODELプロシジャや、SAS/IMLを使用して、多次元正規分布に従う乱数列を生成することも可能です。
詳細については、以下のFAQをご参照ください。

• SAS/ETSで多次元正規分布に従う乱数列を生成する方法
• SAS/IMLで多次元正規分布に従う乱数列を生成する方法