エクセルで多変量正規乱数を発生する方法を考えてみたいと思います。
多変量正規乱数
M次元の多変量乱数はM次元空間の点です。その座標はベクトルであり、M個の成分はランダムで、一般的に独立同一に分布ししているわけではありません。多変量正規変数の特別な場合として多変量ガウス密度関数、
N(x|\mu \Sigma)=\frac{1}{(2\pi)^{M/2}}\exp[-\frac{1}{2}(x-\mu)-\Sigma^{-1}(x-\mu)]
があります。$\mu$は分布の平均をあらわすベクトルです。$\Sigma$は対称な正定値行列で、分布の共分散を表します。$\mu$と$\Sigma$で規定される乱数$x$のベクトルを得る方法は平均ゼロ、分散１の独立な乱数ベクトル$y$を生成することからはじめます。$\Sigma$が正定値行列である場合には、それは下三角行列$L$とその転置行列に分解することができます。
$$ \Sigma=LL^T$$
これをコレスキー分解といいます。つぎに、乱数$x$は分散が１の独立した乱数からつぎの様に生成します。
$$ x = Ly+\mu$$
$x_i$は$y_i$の線形結合として作られます。$i$は要素を表します。平均と共分散によって定義された正規乱数の線形結合は正規分布にしたがいます。
M=2 エクセルシート
まずM=2の場合を考えてみます。$x1$と$x2$の相関係数は
$$\rho=\frac{\text{Cov}(x_1,x_2)}{\sqrt{\text{Var}(x_1)}\sqrt{\text{var}(x_2)}}$$
Var, Covは分散と共分散を表します。平均ゼロ、分散１の独立な乱数$y_1$と$y_2$を生成し、$y_1$を$x_1$とし、$x_2=a_1y_1+a_2y_2$とします。
そうすると、
\begin{eqnarray}
\text{Cov}(x_1,x_2)&=& \text{Cov}(y_1,a_1y_1+a_2y_2)\\
&=& a_1\text{Cov}(y_1,y_1)+a_2\text{Cov}(y_1,y_2)\\
&=& a_1
\end{eqnarray}
また、
\begin{eqnarray}
\text{Var}(x_2)&=& \text{Var}(a_1y_1+a_2y_2)\\
&=& a_1^2\text{Var}(y_1)+a_2^2\text{Var}(y_2)\\
&=& a_1^2+a_2^2=1
\end{eqnarray}
よって、
\begin{eqnarray}
\rho&=&\frac{\text{Cov}(x_1,x_2)}{\sqrt{\text{Var}(x_1)}\sqrt{\text{var}(x_2)}}\\
&=&\frac{a_1}{a_1^2+a_2^2}\\
&=&a_1
\end{eqnarray}
この式を$a_2$について解くと
$$a_2=\frac{\sqrt{1-\rho^2}}{\rho}a_1$$
したがって、
$$x_2=a_1y_1+a_2y_2=a_1y1+\frac{\sqrt{1-\rho^2}}{\rho}a_1=\rho y_1+\sqrt{1-\rho^2}y_2$$
となります。
例 $\rho=0.5$のときb=0.87でエクセルで22000個乱数生成した結果相関は0.496
M＝3 エクセルシート
M＝3では、
\begin{eqnarray}
x_1&=&y_1\\
x_2&=&a_{21}y_1+a_{22}y_2\\
x_3&=&a_{31}y_1+a_{32}y_2+a_{33}y3
\end{eqnarray}
とします。相関は３つあります。x1とx2、x2とx3、そしてx1とx3の間の相関です。それぞれを$\rho_{12},\rho_{23},\rho_{13}$
とします。
したがって
先ほどの計算から
$$x_2=\rho_{12} y_1+\sqrt{1-\rho_{12}^2}y_2$$
となります。また、
\begin{eqnarray}
\text{Cov}(x_2,x_3)&=& \text{Cov}(a_{21}y_1+a_{22}y_2,a_{31}y_1+a_{32}y_2+a_{33}y_3)\\
&=& a_{21}a_{31}\text{Cov}(y_1,y_1)+a_{31}a_{32}\text{Cov}(y_1,y_2)+a_{22}a_{32}\text{Cov}(y_2,y_2)\\
&+& a_{21}a_{33}\text{Cov}(y_1,y_3)+a_{22}a_{33}\text{Cov}(y_2,y_3)\\
&=& a_{21}a_{31}+a_{22}a_{32}\\
&=& \rho_{12}\rho_{13}+\sqrt{1-\rho_{12}}a_{32}\\
&=&\rho_{23}
\end{eqnarray}
よって
$$a_{32}=\frac{\rho_{23}-\rho_{12}\rho_{13}}{\sqrt{1-\rho_{12}}}=$$
\begin{eqnarray}
\text{Var}(x_3)&=& \text{Var}(a_{31}y_1+a_{32}y_2+a_{33}y_3)\\
&=& a_{31}^2\text{Var}(y_1)+a_{32}^2\text{Var}(y_2)+a_{33}^2\text{Var}(y_2)\\
&=& a_{31}^2+a_{32}^2+a_{33}^2\\
&=&\rho_{31}^2+\frac{(\rho_{23}-\rho_{12}\rho_{13})^2}{1-\rho_{12}}+a_{33}^2\\
&=&1
\end{eqnarray}
したがって、
$$a_{33}=\sqrt{1-\rho_{31}^2-\frac{(\rho_{23}-\rho_{12}\rho_{13})^2}{1-\rho_{12}}}$$
例 $\rho_{12}=0.5$, $\rho_{13}=-0.1$, $\rho_{23}=0.7$のとき
M=3 コレスキー分解をエクセルVBAで実行
エクセルVBAでコレスキー分解を行います。エクセルシートではなく、エクセルVBAを用います。Numerical Recipes Third edition のCのコードをエクセルVBAに変換したものを使います。詳しくはp.100 2.9 Cholesky Decompositionを参照してください。
Option Explicit
Sub choldc(el() As Double, n As Long)
Dim i As Long, j As Long, k As Long
Dim sum As Double
For i = 1 To n
For j = i To n
sum = el(i, j)
For k = i - 1 To 1 Step -1
sum = sum - el(i, k) * el(j, k)
Next k
If (i = j) Then
If (sum <= 0#) Then
MsgBox "choldc failed"
End If
el(i, i) = Sqr(sum)
Else
el(j, i) = sum / el(i, i)
End If
Next j
Next i
For i = 1 To n
For j = 0 To i - 1
el(j, i) = 0
Next j
Next i
End Sub
Private Sub ボタン2_Click()
Dim n As Long
Dim el() As Double, b() As Double
Dim i As Long, j As Long
n = 3
ReDim el(n, n)
ReDim b(n, n)
For i = 1 To n
For j = 1 To i
If i = j Then
el(j, i) = 1
Else
el(j, i) = Worksheets("sheet3").Cells(i * 2 - 2, j)
End If
Next j
Next i
Call choldc(el, n)
For j = 1 To n
For i = j To n
If i = 1 And j = 1 Then
Else
Worksheets("sheet3").Cells(i * 2 - 2 + 5, j) = el(i, j)
End If
Next i
Next j
エクセルシートのボタンを押すと実行されます。結果はSheet 2と同じです。
参考
多変量正規分布に従う乱数生成の背景にはCholesky分解がいた
Numerical Recipes in CのCholesky分解ルーチンをExcel VBAに移植
Numerical recipes in C