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

[質問]

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

[回答]

SAS 8以降では、VNORMAL Callを利用できます。
また、ROOT関数を用いたCholesky分解に基づいて、行列演算によって生成することもできます。

以下の例ではいずれも、3次元正規分布に従う乱数列を100オブザベーション作成しています。

■ VNORMAL Call

SAS 8以降では、SAS/IMLのVNORMAL Callを使用して生成できます。
なお、これによって生成される乱数列は旧来の乗算型合同法に基づくものです。

                                            /* VNORMAL Call の例 */
  PROC IML;
    nobs=100;                      /* 生成するオブザベーション数 */
    seed=12345;                              /* 乱数系列のシード */
    mean={20 30 40};                               /* 平均の定義 */

    cov={1    1.2  2.25,
         1.2  4    3.3,
         2.25 3.3  9};                       /* 共分散行列の定義 */

                   /* 戻り値の乱数列、平均ベクトル、共分散行列、 */
                    /*  オブザベーション数、シードの順番で指定。 */
    CALL VNORMAL(rv,mean,cov,nobs,seed);

                                    /* SASデータセットへの出力。 */
                       /* データセット mnormal3 が作成されます。 */
    CREATE mnormal3 FROM rv;
    append from rv;
  quit;

■ Cholesky 分解に基づく処理を行なう

Cholesky分解を行なうROOT関数を使用して、次のようなプログラムで生成できます。
なお、SAS9以降でサポートされているRANDSEED CallとRANDGEN Callを使用すると、MTに基づく乱数列を生成できます。

                                         /* SAS9以降のMTによる例 */
  PROC IML;
    nobs=100;
    seed=12345;

    mean={20 30 40};

    cov={1    1.2  2.25,
         1.2  4    3.3,
         2.25 3.3  9};


    CALL RANDSEED(seed);                       /* シードの初期化 */

                             /* 100×3の行列を事前に作成しておく */
    rvn=J(nobs,NCOL(cov),.);

    CALL RANDGEN(rvn,'NORMAL');                /* 正規乱数の生成 */

                           /* Choleskey 分解を関数 ROOT で行なう */
    rv=mean#J(nobs,NCOL(cov),1)+rvn*ROOT(cov);

                                    /* SASデータセットへの出力。 */
                       /* データセット mnormal4 が作成されます。 */
    CREATE mnormal4 FROM rv;
    APPEND FROM rv;
  quit;


                                             /* 乗算型合同法の例 */
  PROC IML;
    nobs=100;
    seed=12345;

    mean={20 30 40};

    cov={1    1.2  2.25,
         1.2  4    3.3,
         2.25 3.3  9};

    rv=mean#J(nobs,NCOL(cov),1)
         + RANNOR(J(nobs,NCOL(cov),seed))*ROOT(cov);

                                    /* SASデータセットへの出力。 */
                       /* データセット mnormal5 が作成されます。 */
    CREATE mnormal5 FROM rv;
    APPEND FROM rv;
  quit;

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

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