生成一定相关性的二元正态分布

摘要

本文讲简要叙述如何生成具有一定相关性 (correlation) 的服从二元正态分布的随机变量。

二元正态分布

二元正态分布概率密度函数

我们知道单变量标准正态分布$X\sim N(0,1)$ 的概率密度函数 (probability density function, pdf) 为
$$f(x)=\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2},-\infty <x<\infty$$。
那么二元正态分布是如何定义的呢？其定义与单变量的正态分布类似，不过引入了两个变量之间的相关性。具体的定义如下：
如果二元变量 $(X,Y)$ 服从二元正态分布，那么 $(X,Y)$ 的联合概率密度分布函数（joint density distribution function）为：
$$\begin{array}{rl}f(x,y)& =\frac{1}{2\pi {\sigma }_X{\sigma }_Y\sqrt{1-{\rho }^2}}\times \\ & \exp \left(-\frac{1}{2(1-{\rho }^2)}\left(\right(\frac{x-{\mu }_X}{{\sigma }_X}{)}^2-2\rho (\frac{x-{\mu }_X}{{\sigma }_X})(\frac{y-{\mu }_Y}{{\sigma }_Y})+(\frac{y-{\mu }_Y}{{\sigma }_Y}{)}^2\right)\end{array}$$
在公式中，我们称 ${\mu }_X,{\mu }_Y$ 是二元正态分布的均值（期望），${\sigma }_X^2,{\sigma }_Y^2$ 是二元正态分布的方差，$\rho$ 为二元正态分布的相关性。

二元正态分布随机数的生成

现在我们想生成两个服从标准正态分布的随机数，并且要求这两个随机数的相关性是一个给定的数 $\rho$。我们应该如何做呢？

如果我们想用一个程序生成服从标准正态分布的随机数，我们可以用Python numpy包中的 np.random.normal(0, 1, n)，来生成 $n$ 个服从 $N(0,1)$ 的随机数。但是想要生成具有一定相关性的两个服从标准正态分布的随机数，我们就不能用简单得用两次np.random.normal(0, 1, n)，因为这样我们只是生成了两个独立的随机变量，它们的相关性为 0。为了生成具有一定相关性的二元正态分布，这里我们采用构造线性组合的方法。具体如下。

假设 $A,B$ 是两个独立的标准正态分布。 令
$$X=\alpha A+\beta B,$$ $$Y=\gamma A+\delta B$$。
$(\alpha ,\beta ,\gamma ,\delta )$ 是四个待确定的参数。我们希望找到这样的$(\alpha ,\beta ,\gamma ,\delta )$, 使得 $X\sim N(0,1)$， $Y\sim N(0,1)$，并且$\mathrm{C}\mathrm{o}\mathrm{r}(\mathrm{X},\mathrm{Y})=\rho$。我们知道两个独立的正态分布相加依然是一个正态分布。我们有 $X\sim N(0,{\alpha }^2+{\beta }^2)$， $Y\sim N(0,{\gamma }^2+{\delta }^2)$。所以，我们须要${\alpha }^2+{\beta }^2=1$， ${\gamma }^2+{\delta }^2=1$。为了简化运算，我们设 $\alpha =\cos \theta ,\beta =\sin \theta$；$\gamma =\sin \theta ,\delta =\cos \theta$。这样我们的$X,Y$ 就分别服从标准正态分布。下面我们只需要找到 $\theta$，使得$\mathrm{C}\mathrm{o}\mathrm{r}(\mathrm{X},\mathrm{Y})=\rho$。

接下来我们来看条件 $\mathrm{C}\mathrm{o}\mathrm{r}(\mathrm{X},\mathrm{Y})=\rho$。$\mathrm{C}\mathrm{o}\mathrm{r}(\mathrm{X},\mathrm{Y})=\frac{\text{Cov}(\mathrm{X},\mathrm{Y})}{\sqrt{\text{Var}(\mathrm{X})\text{Var}(\mathrm{Y})}}$。因为 $X=(\cos \theta )A+(\sin \theta )B\sim N(0,1)$，所以 $\text{Var}(X)=1$。同样的，$\text{Var}(Y)=1$。

所以我们只须要 $\text{Cov}(X,Y)=\rho$。

而
$$\begin{array}{rl}\text{Cov}(X,Y)& =\text{Cov}((\cos \theta )A+(\sin \theta )B,(\sin \theta )A+(\cos \theta )B)\\ & =\cos \theta \sin \theta +\sin \theta \cos \theta \\ & =\sin 2\theta \end{array}$$。

注意，在上面的计算中我们用到了 $\mathrm{C}\mathrm{o}\mathrm{r}(\mathrm{A},\mathrm{B})=0$。

所以我们有 $\sin 2\theta =\rho$。我们选取 $\theta =\frac{\arcsin \rho }{2}$。这样，我们通过
$$\{\begin{array}{l}X=(\cos \theta )A+(\sin \theta )B\\ Y=(\sin \theta )A+(\cos \theta )B\end{array}$$
计算得到的 $X,Y$ 就是相关性为 $\rho$ 的标准正态分布。

程序实现

我们用Python 来实现上述方法。

import numpy as np
import matplotlib.pyplot as plt

class bivariateNormal:
    
    def __init__(self, rho: 'float', m: int):
        """
        Suppose we want to generate a pair of 
        random variables X, Y, with X ~ N(0, 1), 
        Y ~ N(0, 1), and Cor(X, Y) = rho. m is 
        the number of data pairs we want to generate.
        """
        self.rho = rho
        self.m = m
    
    def generateBivariate(self) -> 'tuple(np.array, np.array)':
        """
        Generate two random variables X, Y, with X ~ N(0, 1), 
        Y ~ N(0, 1), and Cor(X, Y) = rho. 
        self.m is the number of sample points we generated.
        We return a tuple (X, Y). 
        """
        theta = np.arcsin(self.rho) / 2
        A = np.random.normal(0, 1, self.m)
        B = np.random.normal(0, 1, self.m)
        X = np.cos(theta) * A + np.sin(theta) * B
        Y = np.sin(theta) * A + np.cos(theta) * B
        return X, Y

我们试着生成m = 1000个(X, Y)，如下：

m = 10 ** 3
rho = -0.4
a = bivariateNormal(rho, m)
X, Y = a.generateBivariate()
np.corrcoef(X, Y)

得到的结果为：

array([[ 1.        , -0.41408183],
      [-0.41408183,  1.        ]])

即X 与 Y 的相关性为-0.414，与其理论值-0.4 很接近。

多元正态分布的情况

生成服从$N(\mu ,\Sigma )$的$n$ 元正态分布

假设我们想要生成服从$N(\mu ,\mathbf{\Sigma })$ 的 $n$ 元正态分布 $\mathbf{X}=(x_1,x_2,\cdots ,x_n{)}^T$ （$\mathbf{X}$ 是一个列向量，其每个分量 $x_i$ 均为一个随机变量）。我们称 $\mu$ 是 $\mathbf{X}$ 的期望，$\Sigma$ 是 $X$ 的协方差矩阵。

想要生成这样一个随机变量向量， 我们可以采取如下步骤 [1]：

1. 生成 $n$ 个独立的标准正态分布$z_1,z_2,\cdots ,z_n$，记为$\mathbf{Z}$。
2. 找到矩阵 $C$，满足 $C{C}^T=\mathbf{\Sigma }$。

于是 $\mathbf{X}=\mu +C\mathbf{Z}$ 即为满足要求的 $n$ 元随机数。

具体证明的过程可以参考下面的定理：

假设 $\mathbf{X}$ 是一个 $n$-dimensional 的随机变量向量，$\mathbf{X}$ 的期望向量是 $\mu$，$\mathbf{X}$ 的协方差矩阵为 $\Sigma$ 。假设 $\mathbf{C}$ 是一个 $k\times n$ 的实矩阵，$\mathbf{b}$ 是一个 $k\times 1$ 的列向量。则 $\mathbf{Y}=\mathbf{C}\mathbf{X}+\mathbf{b}$ 的期望向量是 $\mathbf{C}\mu +\mathbf{b}$，协方差矩阵为 ${\Sigma }_{\mathbf{Y}}=\mathbf{C}\Sigma {\mathbf{C}}^T$。

多元情况的程序实现

from scipy.linalg import cholesky

class generate_correlated_normal:
    
    def __init__(self, n: int, mu: 'np.array', Sigma: 'np.ndarray', m: int):
        """
        n is the dimension of the random vector.
        mu is the expectation vector of the random vector. 
        mu.shape = (n, 1). 
        Sigma is the covariance matrix of the random sample
        vector we want to generate. m is the number of 
        vectors we want to generate. 
        """
        self.mu = mu.reshape(n, 1)
        self.Sigma = Sigma
        self.m = m
    
    def generate_random_vectors(self) -> 'np.ndarray':
        """
        Generate m random vectors with expectation mu and 
        covariance matrix Sigma.
        The returned value is of type np.ndarray, of shape
        (n, m). n is the dimension of the random vector. 
        """
        C = cholesky(self.Sigma, lower=True)
        n = self.Sigma.shape[0] # n is the dimension of the random vector
        Z = np.random.normal(0, 1, (n, self.m))
        return self.mu + np.dot(C, Z)

n = 3
mu = np.array([0, 0, 0])
Sigma = np.array([[1, 0.2, 0.3], [0.2, 1, 0.4], [0.3, 0.4, 1]])
m = 10 ** 4
b = generate_correlated_normal(n, mu, Sigma, m)
X = b.generate_random_vectors()

np.corrcoef(X[0,:], X[1,:])

array([[1. , 0.2136722],
[0.2136722, 1. ]])

np.corrcoef(X[0,:], X[2,:])

array([[1. , 0.31754401],
[0.31754401, 1. ]])

np.corrcoef(X[1,:], X[2,:])

array([[1. , 0.39934556],
[0.39934556, 1. ]])

可以看出，生成的随机向量元素之间的相关性是符合预期的 [2]。

参考文献

[1] Math stackexchange, Generate Correlated Normal Random Variables

[2] Correlated Random Samples