[PYTHON] Random number generation according to multivariate normal distribution using Cholesky decomposition of covariance matrix

Table of contents

1. Advent Calendar 2020/12/10
2. Multivariate normal distribution
3. Random numbers following a bivariate normal distribution
4. Random numbers following a multivariate normal distribution

Advent Calendar 2020/12/10

In the class of securities investment: theory and practice, which was offered at the Graduate School of Economics of the University of Tokyo (common to all faculties), I tried to write about random number generation according to the multivariate normal distribution using the Cholesky decomposition of the covariance matrix, which was a report subject. think. Wikipedia also has a light explanation, but I wasn't sure what kind of idea this method was devised in the first place, so I'm trying to explain it myself.

Multivariate normal distribution

When the $ n $ dimensional variable $ \ boldsymbol {x} $ follows a multivariate normal distribution with mean $ \ boldsymbol {\ mu} $ and covariance matrix $ \ boldsymbol {\ Sigma} $, the distribution function $ f (\) boldsymbol {x}) $ is $$ \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} $$

Random numbers following a bivariate normal distribution

Let's start with the bivariate (two-dimensional) $ n = 2 $. Let the mean of $ i (i = x, y) $ be $ \ mu_i $, the variance be $ \ sigma_i ^ 2 $, and the correlation coefficient be $ \ rho $. The bivariate normal distribution function that $ x, y $ follows is

    \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}

It will be. So, after generating $ x $ to follow the normal distribution $ \ mathrm {N} (\ mu_x, \ sigma_x ^ 2) $, $ y $ follows the normal distribution $ \ mathrm {N} (\ mu_y + \ frac {\ sigma_y} {\ sigma_x} \ rho (x- \ mu_x), (1- \ rho ^ 2) \ sigma_y ^ 2) $ You can see that it should be generated in. As a program example

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()

Random numbers following a multivariate normal distribution

Even if the above method is bivariate, the formula becomes complicated as the dimension increases, and it is difficult to program for random number generation. Also, it cannot respond to changes in dimensions. So, let's take a look at the distribution function to see if we can generate random numbers in a more comfortable and flexible way. Using $ n $ dimensional vector $ z $ and $ n \ times n $ matrix $ \ boldsymbol {A} $ $$ \begin{equation} \boldsymbol{x}-\boldsymbol{\mu}=\boldsymbol{A}\boldsymbol{z} \end{equation} $$ And let the distribution function of $ \ boldsymbol {z} $ be $ g (\ boldsymbol {z}) $. From the transformation of the probability density function accompanying the change of variables

    \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}

It will be. if,

    \begin{equation}
        \boldsymbol{A}^\mathrm{T}\boldsymbol{\Sigma}^{-1}\boldsymbol{A}=\boldsymbol{1} 
    \end{equation}

If so, $ \begin{equation} |\boldsymbol{A}|=|\boldsymbol{\Sigma}|^\frac{1}{2} \end{equation} $ So $ \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} $ And $ \ boldsymbol {z} $ shows that each element follows an independent standard normal distribution. for that purpose

    \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}

Must be expressed as. If the covariance matrix is ​​a definite-value symmetric matrix, then such $ \ boldsymbol {A} $ exists and its value can be determined by Cholesky decomposition of the covariance matrix. That's why Cholesky decomposition is used. The great thing about this method is that it's fairly easy to program and it's easy to adapt to changing dimensions. The following is a trivariate program.

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()