目次
Advent Calendar 2020/12/10
東大経済学研究科開講（学部共通）の証券投資：理論と実践という授業にて、レポート課題になっていた、共分散行列のコレスキー分解を用いた多変量正規分布に従う乱数生成について書こうと思います。Wikipediaにも軽く説明はありますが、そもそもどういう発想でこの方法が考案されたのか不明だったので自分なりに説明しようという試みです。
多変量正規分布
$ n $次元の変数$ \boldsymbol{x} $が、平均$ \boldsymbol{\mu} $、共分散行列$ \boldsymbol{\Sigma} $の多変量正規分布に従うとき、分布関数$ f(\boldsymbol{x}) $は
$$
\begin{equation}
f(\boldsymbol{x})=\frac{1}{(2\pi)^\frac{n}{2}|\boldsymbol{\Sigma}|^\frac{1}{2}}
\exp\left[-\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^\mathrm{T}\boldsymbol{\Sigma}^{-1}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)\right]
\end{equation}
$$
二変量正規分布に従う乱数
まずは二変量(二次元)$ n=2 $から見てみます。
$ i(i=x,y) $の平均を$ \mu_i $、分散を$ \sigma_i^2 $、相関係数を$ \rho $とすると、
$ x,y $が従う二変量正規分布関数は
\begin{align}
f(x,y)&=\frac{1}{2\pi\sigma_x\sigma_y\sqrt{1-\rho^2}}
\exp\left[-\frac{1}{2(1-\rho^2)}
\left(\frac{(x-\mu_x)^2}{\sigma_x^2}
-2\rho\frac{(x-\mu_x)(y-\mu_y)}{\sigma_x\sigma_y}
+\frac{(y-\mu_y)^2}{\sigma_y^2}\right)\right] \\
&=
\frac{1}{2\pi\sigma_x\sigma_y\sqrt{1-\rho^2}}
\exp\left[-\frac{(x-\mu_x)^2}{2\sigma_x^2}
-\frac{(y-(\mu_y+\frac{\sigma_y}{\sigma_x}\rho(x-\mu_x)))^2}{2(1-\rho^2)\sigma_y^2}\right]
\end{align}
となります。
なので、$ x $を正規分布$ \mathrm{N}(\mu_x,\sigma_x^2) $に従うように発生させたのち、
$ y $は正規分布$ \mathrm{N}(\mu_y+\frac{\sigma_y}{\sigma_x}\rho(x-\mu_x),(1-\rho^2)\sigma_y^2) $に従うように発生させればよいことがわかります。
プログラム例として
binorm.py
import numpy as np
import matplotlib.pyplot as plt
n = 100
mux, muy = 40, 60
sigmax, sigmay = 10, 5
rho = 0.7
randomx = np.random.normal(mux, sigmax, [n])
muy2 = muy + sigmay/sigmax * rho * (randomx - mux)
sigmay2 = ((1 - rho**2) * sigmay*2)**0.5
randomy = np.random.normal(muy2, sigmay2)
plt.figure()
plt.scatter(randomx,randomy)
plt.xticks([mux-2*sigmax, mux-sigmax, mux, mux+sigmax, mux+2*sigmax]
,[r'$\mu_x-2\sigma_x$', r'$\mu_x-\sigma_x$', r'$\mu_x$', r'$\mu_x+\sigma_x$', r'$\mu_x+2\sigma_x$'])
plt.xlabel('x')
plt.yticks([muy-2*sigmay, muy-sigmay, muy, muy+sigmay, muy+2*sigmay]
,[r'$\mu_y-2\sigma_y$', r'$\mu_y-\sigma_y$', r'$\mu_y$', r'$\mu_y+\sigma_y$', r'$\mu_y+2\sigma_y$'])
plt.ylabel('y')
plt.savefig('binorm.png')
plt.close()
多変量正規分布に従う乱数
上述の方法は二変量ならまだしも、次元が大きくなるにつれて式が複雑になり、乱数発生のためのプログラムを組むのに一苦労です。また、次元の変化に対応できません。そこで、もっと楽に柔軟な方法で乱数を発生させられないか、分布関数を眺めてみます。$ n $次元ベクトル$ z $と$ n \times n $行列$ \boldsymbol{A} $を用いて
$$
\begin{equation}
\boldsymbol{x}-\boldsymbol{\mu}=\boldsymbol{A}\boldsymbol{z}
\end{equation}
$$
と置き、$ \boldsymbol{z} $の分布関数を$ g(\boldsymbol{z}) $とします。変数変換に伴う確率密度関数の変換より
\begin{align}
g(\boldsymbol{z})d\boldsymbol{z}&=f(\boldsymbol{x})d\boldsymbol{x} \\
g(\boldsymbol{z})&=f(\boldsymbol{x})\left|\frac{d\boldsymbol{x}}{d\boldsymbol{z}}\right| \\
&=f(\boldsymbol{A}\boldsymbol{z}+\boldsymbol{\mu})|\boldsymbol{A}| \\
&=\frac{|\boldsymbol{A}|}{(2\pi)^\frac{n}{2}|\boldsymbol{\Sigma}|^\frac{1}{2}}
\exp\left[-\frac{1}{2}\left(\boldsymbol{A}\boldsymbol{z}\right)^\mathrm{T}
\boldsymbol{\Sigma}^{-1}\left(\boldsymbol{A}\boldsymbol{z}\right)\right] \\
&=\frac{|\boldsymbol{A}|}{(2\pi)^\frac{n}{2}|\boldsymbol{\Sigma}|^\frac{1}{2}}
\exp\left[-\frac{1}{2}\boldsymbol{z}^\mathrm{T}\boldsymbol{A}^\mathrm{T}
\boldsymbol{\Sigma}^{-1}\boldsymbol{A}\boldsymbol{z}\right]
\end{align}
となります。もし、
\begin{equation}
\boldsymbol{A}^\mathrm{T}\boldsymbol{\Sigma}^{-1}\boldsymbol{A}=\boldsymbol{1}
\end{equation}
となれば、
$$
\begin{equation}
|\boldsymbol{A}|=|\boldsymbol{\Sigma}|^\frac{1}{2}
\end{equation}
$$
なので、
$$
\begin{equation}
g(\boldsymbol{z})=\frac{1}{(2\pi)^\frac{n}{2}}
\exp\left[-\frac{1}{2}\boldsymbol{z}^\mathrm{T}\boldsymbol{z}\right]
\end{equation}
$$
となり、$ \boldsymbol{z} $は各要素が独立な標準正規分布に従っていることがわかります。そのためには
\begin{eqnarray}
\boldsymbol{\Sigma}^{-1}&=&\left(\boldsymbol{A}^\mathrm{T}\right)^{-1}\boldsymbol{A}^{-1} \\
\boldsymbol{\Sigma}&=&\left[\left(\boldsymbol{A}^\mathrm{T}\right)^{-1}\boldsymbol{A}^{-1}\right]^{-1} \\
&=&\boldsymbol{A}\boldsymbol{A}^\mathrm{T}
\end{eqnarray}
と表されることが必要です。共分散行列が正定値対称行列である場合には、このような$ \boldsymbol{A} $は存在し、その値は共分散行列をコレスキー分解することによって求めることができます。以上が、コレスキー分解が用いられる理由です。この方法の素晴らしいところはプログラムがかなり簡単にかけ、次元の変化にも対応しやすい点でしょうか。以下は三変量のプログラムです。
multinorm.py
import numpy as np
import matplotlib.pyplot as plt
n = 300
dim = 3
mux, muy, muz = 40, 60, 80
mu = np.array([[mux]
,[muy]
,[muz]])
sigmax, sigmay, sigmaz = 10, 5, 15
sigmaxy, sigmayz, sigmazx = 0.3 * sigmax * sigmay, 0.5 * sigmay * sigmaz, 0.7 * sigmaz * sigmax
Cov = np.array([[sigmax**2, sigmaxy, sigmazx]
,[sigmaxy, sigmay**2, sigmayz]
,[sigmazx, sigmayz, sigmaz**2]])
upTriang = np.linalg.cholesky(Cov)
randomz = np.random.randn(dim, n)
randomx, randomy, randomz = mu + np.dot(upTriang, randomz)
fig = plt.figure(figsize=(10,30), tight_layout=True)
ax1, ax2, ax3 = fig.add_subplot(311),fig.add_subplot(312), fig.add_subplot(313)
ax1.scatter(randomx,randomy)
ax1.set_xticks([mux-2*sigmax, mux-sigmax, mux, mux+sigmax, mux+2*sigmax])
ax1.set_xticklabels([r'$\mu_x-2\sigma_x$', r'$\mu_x-\sigma_x$', r'$\mu_x$', r'$\mu_x+\sigma_x$', r'$\mu_x+2\sigma_x$'])
ax1.set_xlabel('x')
ax1.set_yticks([muy-2*sigmay, muy-sigmay, muy, muy+sigmay, muy+2*sigmay])
ax1.set_yticklabels([r'$\mu_y-2\sigma_y$', r'$\mu_y-\sigma_y$', r'$\mu_y$', r'$\mu_y+\sigma_y$', r'$\mu_y+2\sigma_y$'])
ax1.set_ylabel('y')
ax1.grid()
ax2.scatter(randomy,randomz)
ax2.set_xticks([muy-2*sigmay, muy-sigmay, muy, muy+sigmay, muy+2*sigmay])
ax2.set_xticklabels([r'$\mu_y-2\sigma_y$', r'$\mu_y-\sigma_y$', r'$\mu_y$', r'$\mu_y+\sigma_y$', r'$\mu_y+2\sigma_y$'])
ax2.set_xlabel('y')
ax2.set_yticks([muz-2*sigmaz, muz-sigmaz, muz, muz+sigmaz, muz+2*sigmaz])
ax2.set_yticklabels([r'$\mu_z-2\sigma_z$', r'$\mu_z-\sigma_z$', r'$\mu_z$', r'$\mu_z+\sigma_z$', r'$\mu_z+2\sigma_z$'])
ax2.set_ylabel('z')
ax2.grid()
ax3.scatter(randomz,randomx)
ax3.set_xticks([muz-2*sigmaz, muz-sigmaz, muz, muz+sigmaz, muz+2*sigmaz])
ax3.set_xticklabels([r'$\mu_z-2\sigma_z$', r'$\mu_z-\sigma_z$', r'$\mu_z$', r'$\mu_z+\sigma_z$', r'$\mu_z+2\sigma_z$'])
ax3.set_xlabel('z')
ax3.set_yticks([mux-2*sigmax, mux-sigmax, mux, mux+sigmax, mux+2*sigmax])
ax3.set_yticklabels([r'$\mu_x-2\sigma_x$', r'$\mu_x-\sigma_x$', r'$\mu_x$', r'$\mu_x+\sigma_x$', r'$\mu_x+2\sigma_x$'])
ax3.set_ylabel('x')
ax3.grid()
fig.savefig('multinorm.png')
plt.close()