二変量正規乱数発生のアルゴリズム

前回の記事から、二変量の正規乱数の発生方法が欲しいなと思いましたので、エッセンスだけ取りまとめました。


ボックスミュラーのアルゴリズム - nishiru3の日記

検索をかけたら、非常に良い資料が出てきました。
下記URLのうち、アルゴリズムに関係ある部分だけをピックアップしております。
http://www012.upp.so-net.ne.jp/doi/sas/simulation/multi_co/bivariate_normal.pdf


二変量の場合、平均値と分散だけではなく、共分散も必要になります。xとyの関係も必要ということです。

\begin{eqnarray}
\left(
\begin{array}{cc}
      x \\
      y 
\end{array}
\right)\sim N
\left(
\begin{array}{cc}
     \left(
     \begin{array}{cc}
      \mu_x \\
      \mu_y
	\end{array}
     \right)
     ,
     \left(
     \begin{array}{cc}
     {\sigma_x}^2 & \sigma_{yx} \\
	\sigma_{xy} & {\sigma_y}^2  
	\end{array}
     \right)
\end{array}
\right)
\end{eqnarray}

まず下記の条件でxの正規乱数を発生させます。
\begin{eqnarray}
x \sim N\left(\mu_x,{\sigma_x}^2\right)
\end{eqnarray}
ここで、xの平均を\mu_x、分散を{\sigma_x}^2とします。
これは、xに関する正規乱数になります。発生方法はボックスミューラー法を用いることにします。
さて、yについても同様に正規乱数を発生させる必要がありますが、そのままだと、yも独立に正規乱数を発生させることになります。

そこで、yの発生方法は以下のように考えます。

\begin{eqnarray}
y \sim N\left(\mu_y+\rho \frac{\sigma_y}{\sigma_x}(x-\mu_x),(1-\rho^2){\sigma_y}\right)
\end{eqnarray}
ここで、当初{\sigma_{xy}=\rho \sigma_x \sigma_y}と考えます。\rhoはx、yの相関係数です。

ちなみに、x,yの平均値を0、相関係数を\rho=0とすると、標準正規分布となるようです。

多変量であっても、各変量の相関係数を設定してあげれば、同様に発生させることができると思います。ただ、ちょっと煩雑になるかもしれません。