如果我们想用一个程序生成服从标准正态分布的随机数，我们可以用Python numpy包中的 np.random.normal(0, 1, n)，来生成
n
n
n 个服从
N
(
0
,
1
)
N(0, \, 1)
N(0,1) 的随机数。但是想要生成具有一定相关性的两个服从标准正态分布的随机数，我们就不能用简单得用两次np.random.normal(0, 1, n)，因为这样我们只是生成了两个独立的随机变量，它们的相关性为 0。为了生成具有一定相关性的二元正态分布，这里我们采用构造线性组合的方法。具体如下。
假设
A
,
B
A, \, B
A,B 是两个独立的标准正态分布。 令
X
=
α
A
+
β
B
,
X = \alpha A + \beta B,
X=αA+βB,
Y
=
γ
A
+
δ
B
Y = \gamma A + \delta B
Y=γA+δB。
(
α
,
β
,
γ
,
δ
)
(\alpha, \, \beta, \, \gamma, \, \delta)
(α,β,γ,δ) 是四个待确定的参数。我们希望找到这样的
(
α
,
β
,
γ
,
δ
)
(\alpha, \, \beta, \, \gamma, \, \delta)
(α,β,γ,δ), 使得
X
∼
N
(
0
,
1
)
X \sim N(0, \, 1)
X∼N(0,1)，
Y
∼
N
(
0
,
1
)
Y \sim N(0, \, 1)
Y∼N(0,1)，并且
C
o
r
(
X
,
Y
)
=
ρ
\rm{Cor}(X, \, Y) = \rho
Cor(X,Y)=ρ。我们知道两个独立的正态分布相加依然是一个正态分布。我们有
X
∼
N
(
0
,
α
2
+
β
2
)
X \sim N(0, \, \alpha^2 + \beta^2)
X∼N(0,α2+β2)，
Y
∼
N
(
0
,
γ
2
+
δ
2
)
Y \sim N(0, \, \gamma^2 + \delta^2)
Y∼N(0,γ2+δ2)。所以，我们须要
α
2
+
β
2
=
1
\alpha^2 + \beta^2 = 1
α2+β2=1，
γ
2
+
δ
2
=
1
\gamma^2 + \delta^2 = 1
γ2+δ2=1。为了简化运算，我们设
α
=
cos
θ
,
β
=
sin
θ
\alpha = \cos\theta, \, \beta = \sin \theta
α=cosθ,β=sinθ；
γ
=
sin
θ
,
δ
=
cos
θ
\gamma = \sin \theta, \, \delta = \cos \theta
γ=sinθ,δ=cosθ。这样我们的
X
,
Y
X, \, Y
X,Y 就分别服从标准正态分布。下面我们只需要找到
θ
\theta
θ，使得
C
o
r
(
X
,
Y
)
=
ρ
\rm{Cor}(X, \, Y) = \rho
Cor(X,Y)=ρ。
接下来我们来看条件
C
o
r
(
X
,
Y
)
=
ρ
\rm{Cor}(X, \, Y) = \rho
Cor(X,Y)=ρ。
C
o
r
(
X
,
Y
)
=
Cov
(
X
,
Y
)
Var
(
X
)
Var
(
Y
)
\displaystyle \rm{Cor}(X, \, Y) = \frac{\text{Cov}(X, \, Y)}{\sqrt{\text{Var}(X) \text{Var}(Y)}}
Cor(X,Y)=Var(X)Var(Y)Cov(X,Y)。因为
X
=
(
cos
θ
)
A
+
(
sin
θ
)
B
∼
N
(
0
,
1
)
X = (\cos \theta) A + (\sin \theta) B \sim N(0, \, 1)
X=(cosθ)A+(sinθ)B∼N(0,1)，所以
Var
(
X
)
=
1
\text{Var}(X) = 1
Var(X)=1。同样的，
Var
(
Y
)
=
1
\text{Var}(Y) = 1
Var(Y)=1。
所以我们只须要
Cov
(
X
,
Y
)
=
ρ
\text{Cov}(X, \, Y) = \rho
Cov(X,Y)=ρ。
而
Cov
(
X
,
Y
)
=
Cov
(
(
cos
θ
)
A
+
(
sin
θ
)
B
,
(
sin
θ
)
A
+
(
cos
θ
)
B
)
=
cos
θ
sin
θ
+
sin
θ
cos
θ
=
sin
2
θ
\begin{aligned} \text{Cov}(X, \, Y) &= \text{Cov} \Big( (\cos \theta) A + (\sin \theta) B, (\sin \theta) A + (\cos \theta) B \Big) \\ &= \cos \theta \sin \theta + \sin \theta \cos \theta \\ &= \sin2\theta \end{aligned}
Cov(X,Y)=Cov((cosθ)A+(sinθ)B,(sinθ)A+(cosθ)B)=cosθsinθ+sinθcosθ=sin2θ。
注意，在上面的计算中我们用到了
C
o
r
(
A
,
B
)
=
0
\rm{Cor}(A, \, B) = 0
Cor(A,B)=0。
所以我们有
sin
2
θ
=
ρ
\sin2\theta = \rho
sin2θ=ρ。我们选取
θ
=
arcsin
ρ
2
\displaystyle \theta = \frac{ \arcsin \rho}{2}
θ=2arcsinρ。这样，我们通过
{
X
=
(
cos
θ
)
A
+
(
sin
θ
)
B
Y
=
(
sin
θ
)
A
+
(
cos
θ
)
B
\begin{cases} X = (\cos \theta) A + (\sin \theta) B \\ Y = (\sin \theta) A + (\cos \theta) B \end{cases}
{X=(cosθ)A+(sinθ)BY=(sinθ)A+(cosθ)B 计算得到的
X
,
Y
X, \, Y
X,Y 就是相关性为
ρ
\rho
ρ 的标准正态分布。
程序实现
我们用Python 来实现上述方法。
import numpy as np
import matplotlib.pyplot as plt
classbivariateNormal: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
defgenerateBivariate(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
(
μ
,
Σ
)
N(\mathbf{\mu}, \, \Sigma)
N(μ,Σ)的
n
n
n 元正态分布
假设我们想要生成服从
N
(
μ
,
Σ
)
N(\mathbf{\mu}, \, \mathbf{\Sigma})
N(μ,Σ) 的
n
n
n 元正态分布
X
=
(
x
1
,
x
2
,
⋯
,
x
n
)
T
\mathbf{X} = (x_1, \, x_2, \cdots, \, x_n)^T
X=(x1,x2,⋯,xn)T （
X
\mathbf{X}
X 是一个列向量，其每个分量
x
i
x_i
xi 均为一个随机变量）。我们称
μ
\mu
μ 是
X
\mathbf{X}
X 的期望，
Σ
\Sigma
Σ 是
X
X
X 的协方差矩阵。
想要生成这样一个随机变量向量， 我们可以采取如下步骤 [1]：
生成
n
n
n 个独立的标准正态分布
z
1
,
z
2
,
⋯
,
z
n
z_1, \, z_2, \cdots, \, z_n
z1,z2,⋯,zn，记为
Z
\mathbf{Z}
Z。
找到矩阵
C
C
C，满足
C
C
T
=
Σ
CC^T = \mathbf{\Sigma}
CCT=Σ。
于是
X
=
μ
+
C
Z
\mathbf{X} = \mu + C \mathbf{Z}
X=μ+CZ 即为满足要求的
n
n
n 元随机数。
具体证明的过程可以参考下面的定理：
假设
X
\mathbf{X}
X 是一个
n
n
n-dimensional 的随机变量向量，
X
\mathbf{X}
X 的期望向量是
μ
\mu
μ，
X
\mathbf{X}
X 的协方差矩阵为
Σ
\Sigma
Σ 。假设
C
\mathbf{C}
C 是一个
k
×
n
k \times n
k×n 的实矩阵，
b
\mathbf{b}
b 是一个
k
×
1
k \times 1
k×1 的列向量。则
Y
=
C
X
+
b
\mathbf{Y} = \mathbf{C} \mathbf{X} + \mathbf{b}
Y=CX+b 的期望向量是
C
μ
+
b
\mathbf{C} \mu + \mathbf{b}
Cμ+b，协方差矩阵为
Σ
Y
=
C
Σ
C
T
\Sigma_{\mathbf{Y}} = \mathbf{C} \Sigma \mathbf{C}^T
ΣY=CΣCT。
多元情况的程序实现
from scipy.linalg import cholesky
classgenerate_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
defgenerate_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()
