协方差矩阵的模拟及独立性、相关性的判断

本博客主要解决以下3个问题：
(i)利用矩阵分析的相关理论，模拟产生一个协方差矩阵；

(ii)由此协方差矩阵模拟产生一个多维正态随机变量随机数，并分析该组数据间的独立性；

(iii)如果具有相关性，如何将其转化为不相关？

协方差矩阵的模拟及独立性、相关性的判断

• 1 协方差矩阵的模拟
• • 1.1 协方差矩阵的基础知识
• 2 多元正态随机数的产生
• 3 相关性的判断

1 协方差矩阵的模拟

1.1 协方差矩阵的基础知识

（1）多维随机变量的协方差矩阵

对多维随机变量 X = [ X 1 , X 2 , … , X n ] T \mathbf{X}=\left[X_{1}, X_{2}, \ldots, X_{n}\right]^{T} X=[X1​,X2​,…,Xn​]T，我们往往需要计算各维度之间的协方差，这样协方差就组成了一个 n × n n \times n n×n的矩阵，称为协方差矩阵。协方差矩阵是一个对角矩阵，对角线上的元素是各维度上随机变量的方差。 我们定义协方差为 Σ \Sigma Σ, 矩阵内的元素 Σ i j \Sigma_{i j} Σij​为

Σ i j = cov ⁡ ( X i , X j ) = E [ ( X i − E ( X i ) ) ( X j − E ( X j ) ) ] \Sigma_{i j}=\operatorname{cov}\left(X_{i}, X_{j}\right)=E\left[\left(X_{i}-E\left(X_{i}\right)\right)\left(X_{j}-E\left(X_{j}\right)\right)\right] Σij​=cov(Xi​,Xj​)=E[(Xi​−E(Xi​))(Xj​−E(Xj​))]

协方差矩阵为
Σ = E [ ( X − E ( X ) ) ( X − E ( X ) ) T ] = [ cov ⁡ ( X 1 , X 1 ) cov ⁡ ( X 1 , X 2 ) ⋯ cov ⁡ ( X 1 , X n ) cov ⁡ ( X 2 , X 1 ) cov ⁡ ( X 2 , X 2 ) ⋯ cov ⁡ ( X 2 , X n ) ⋮ ⋮ ⋱ ⋯ cov ⁡ ( X n , X 1 ) cov ⁡ ( X n , X 2 ) ⋯ cov ⁡ ( X n , X n ) ] . \begin{aligned} \Sigma &=E\left[(X-E(X))(X-E(X))^{T}\right] \\ &=\left[\begin{array}{cccc} \operatorname{cov}\left(X_{1}, X_{1}\right) & \operatorname{cov}\left(X_{1}, X_{2}\right) & \cdots & \operatorname{cov}\left(X_{1}, X_{n}\right) \\ \operatorname{cov}\left(X_{2}, X_{1}\right) & \operatorname{cov}\left(X_{2}, X_{2}\right) & \cdots & \operatorname{cov}\left(X_{2}, X_{n}\right) \\ \vdots & \vdots & \ddots & \cdots \\ \operatorname{cov}\left(X_{n}, X_{1}\right) & \operatorname{cov}\left(X_{n}, X_{2}\right) & \cdots & \operatorname{cov}\left(X_{n}, X_{n}\right) \end{array}\right] . \end{aligned} Σ​=E[(X−E(X))(X−E(X))T]=⎣⎢⎢⎢⎡​cov(X1​,X1​)cov(X2​,X1​)⋮cov(Xn​,X1​)​cov(X1​,X2​)cov(X2​,X2​)⋮cov(Xn​,X2​)​⋯⋯⋱⋯​cov(X1​,Xn​)cov(X2​,Xn​)⋯cov(Xn​,Xn​)​⎦⎥⎥⎥⎤​.​

（2）样本的协方差矩阵

假设数据集 T = { x i } i = 1 m T=\left\{x_{i}\right\}_{i=1}^{m} T={xi​}i=1m​表示 m m m个样本，每个样本表示为 X i = ( x i 1 , x i 2 , … , x i n ) T {{X}_{i}}={{\left( {{x}_{i1}},{{x}_{i2}},\ldots ,{{x}_{in}} \right)}^{T}} Xi​=(xi1​,xi2​,…,xin​)T。所有样本可以组成一个 m × n m×n m×n的矩阵。

X m × n = [ x 11 x 12 ⋯ x 1 n x 21 x 22 ⋯ x 2 n ⋮ ⋮ ⋮ ⋮ x m 1 x m 2 ⋯ x m n ] = [ c 1 , c 2 , … , c n ] X_{m \times n}=\left[\begin{array}{cccc}x_{11} & x_{12} & \cdots & x_{1 n} \\ x_{21} & x_{22} & \cdots & x_{2 n} \\ \vdots & \vdots & \vdots & \vdots \\ x_{m 1} & x_{m 2} & \cdots & x_{m n}\end{array}\right]=\left[c_{1}, c_{2}, \ldots, c_{n}\right] Xm×n​=⎣⎢⎢⎢⎡​x11​x21​⋮xm1​​x12​x22​⋮xm2​​⋯⋯⋮⋯​x1n​x2n​⋮xmn​​⎦⎥⎥⎥⎤​=[c1​,c2​,…,cn​]

每一行代表一个对象，每一类代表一个维度，协方差矩阵是求维度之间的相关性，而不是对象之间的，所以协方差矩阵的大小与维度相同。 C i C_{i} Ci​表示第 i i i维的随机变量。

假设 x ˉ = ( x ˉ 1 , x ˉ 2 , … , x ˉ n ) \bar{x}=\left(\bar{x}_{1}, \bar{x}_{2}, \ldots, \bar{x}_{n}\right) xˉ=(xˉ1​,xˉ2​,…,xˉn​), 则有 E ( c i ) = x i ‾ E\left( {{c}_{i}} \right)=\overline{{{x}_{i}}} E(ci​)=xi​​。用 Σ \Sigma Σ表示样本的协方差矩阵，则有
Σ = [ cov ⁡ ( c 1 , c 1 ) cov ⁡ ( c 1 , c 2 ) ⋯ cov ⁡ ( c 1 , c n ) cov ⁡ ( c 2 , c 1 ) cov ⁡ ( c 2 , c 2 ) ⋯ cov ⁡ ( c 2 , c n ) ⋮ ⋮ ⋱ ⋯ cov ⁡ ( c n , c 1 ) cov ⁡ ( c n , c 2 ) ⋯ cov ⁡ ( c n , c n ) ] \Sigma=\left[\begin{array}{cccc} \operatorname{cov}\left(c_{1}, c_{1}\right) & \operatorname{cov}\left(c_{1}, c_{2}\right) & \cdots & \operatorname{cov}\left(c_{1}, c_{n}\right) \\ \operatorname{cov}\left(c_{2}, c_{1}\right) & \operatorname{cov}\left(c_{2}, c_{2}\right) & \cdots & \operatorname{cov}\left(c_{2}, c_{n}\right) \\ \vdots & \vdots & \ddots & \cdots \\ \operatorname{cov}\left(c_{n}, c_{1}\right) & \operatorname{cov}\left(c_{n}, c_{2}\right) & \cdots & \operatorname{cov}\left(c_{n}, c_{n}\right) \end{array}\right] Σ=⎣⎢⎢⎢⎡​cov(c1​,c1​)cov(c2​,c1​)⋮cov(cn​,c1​)​cov(c1​,c2​)cov(c2​,c2​)⋮cov(cn​,c2​)​⋯⋯⋱⋯​cov(c1​,cn​)cov(c2​,cn​)⋯cov(cn​,cn​)​⎦⎥⎥⎥⎤​
进一步可化简为
[ ∑ i = 1 m ( x i 1 − x ˉ 1 ) ( x i 1 − x ˉ 1 ) ∑ i = 1 m ( x i 1 − x ˉ 1 ) ( x i 2 − x ˉ 2 ) ⋯ ∑ i = 1 m ( x i 1 − x ˉ 1 ) ( x i n − x ˉ n ) ∑ i = 1 m ( x j 2 − x ˉ 2 ) ( x i 1 − x ˉ 1 ) ∑ i = 1 n ( x i 2 − x ˉ 2 ) ( x i 2 − x ˉ 2 ) ⋯ ∑ i = 1 m ( x i n − x ˉ n ) ( x i n − x ˉ n ) ⋮ ⋮ ⋱ ⋯ ∑ i = 1 m ( x i n − x ˉ 1 ) ( x i n − x ˉ 1 ) ∑ i = 1 m ( x i n − x ˉ n ) ( x i 2 − x ˉ 2 ) ⋯ ∑ i = 1 m ( x i n − x ˉ n ) ( x i n − x ˉ n ) ] = 1 m − 1 ∑ i = 1 m ( x i − x ˉ ) ( x i − x ˉ ) T \begin{array}{l} \left[\begin{array}{cccc} \sum_{i=1}^{m}\left(x_{i 1}-\bar{x}_{1}\right)\left(x_{i 1}-\bar{x}_{1}\right) & \sum_{i=1}^{m}\left(x_{i 1}-\bar{x}_{1}\right)\left(x_{i 2}-\bar{x}_{2}\right) & \cdots & \sum_{i=1}^{m}\left(x_{i 1}-\bar{x}_{1}\right)\left(x_{i n}-\bar{x}_{n}\right) \\ \sum_{i=1}^{m}\left(x_{j 2}-\bar{x}_{2}\right)\left(x_{i 1}-\bar{x}_{1}\right) & \sum_{i=1}^{n}\left(x_{i 2}-\bar{x}_{2}\right)\left(x_{i 2}-\bar{x}_{2}\right) & \cdots & \sum_{i=1}^{m}\left(x_{i n}-\bar{x}_{n}\right)\left(x_{i n}-\bar{x}_{n}\right) \\ \vdots & \vdots & \ddots & \cdots \\ \sum_{i=1}^{m}\left(x_{i n}-\bar{x}_{1}\right)\left(x_{i n}-\bar{x}_{1}\right) & \sum_{i=1}^{m}\left(x_{i n}-\bar{x}_{n}\right)\left(x_{i 2}-\bar{x}_{2}\right) & \cdots & \sum_{i=1}^{m}\left(x_{i n}-\bar{x}_{n}\right)\left(x_{i n}-\bar{x}_{n}\right) \end{array}\right] \\ =\frac{1}{m-1} \sum_{i=1}^{m}\left(x_{i}-\bar{x}\right)\left(x_{i}-\bar{x}\right)^{T} \end{array} ⎣⎢⎢⎢⎡​∑i=1m​(xi1​−xˉ1​)(xi1​−xˉ1​)∑i=1m​(xj2​−xˉ2​)(xi1​−xˉ1​)⋮∑i=1m​(xin​−xˉ1​)(xin​−xˉ1​)​∑i=1m​(xi1​−xˉ1​)(xi2​−xˉ2​)∑i=1n​(xi2​−xˉ2​)(xi2​−xˉ2​)⋮∑i=1m​(xin​−xˉn​)(xi2​−xˉ2​)​⋯⋯⋱⋯​∑i=1m​(xi1​−xˉ1​)(xin​−xˉn​)∑i=1m​(xin​−xˉn​)(xin​−xˉn​)⋯∑i=1m​(xin​−xˉn​)(xin​−xˉn​)​⎦⎥⎥⎥⎤​=m−11​∑i=1m​(xi​−xˉ)(xi​−xˉ)T​
这里分母为 m − 1 m−1 m−1是因为随机变量的数学期望未知，以样本均值代替，*度减一。

（3）模拟产生一个协方差矩阵

方法一计算步骤：

Step1: 先求X的每一列的均值;

Step2: 计算每列减去每列相应的均值；

Step3: 乘积求和。

covdef <- function(x){
  #方法1
  x<-as.matrix(x)
  n<-nrow(x)
  mx<-diag(1,n)-matrix(1,n,n)/n
  cov_matrix<-(t(x)%*%mx%*%x)/(n-1)
  return(cov_matrix)
}

方法二

covfun <- function(x){
  #方法2
  #数据的长度
  m <- nrow(x)
  #变量的个数
  n <- ncol(x)
  #求每一列的平均值
  colmean <- colMeans(x)
  colx <-matrix(nrow =m,ncol = n )
  cov_matrix <-matrix(nrow =n,ncol = n )
  for(i in 1:m){
    for(j in 1:n){
      colx[i,j] <- x[i,j]-colmean[j] 
    }
  }
  cov_matrix <- t(colx)%*%(colx)/(m-1)
  return(cov_matrix)
}

函数的调用

#函数的调用
data <- matrix(c(1,2,3,2,1,3,3,1,2),ncol = 3,nrow  = 3)
#方法一
source('D:/RStudio/R_code/analysis_stastic/covfun.R')
covfun(data)
#方法二
source('D:/RStudio/R_code/analysis_stastic/covdef.R')
covdef(data)

#方法三
#R自带函数包
cov(data)   

结果显示

> covfun(data)
     [,1] [,2] [,3]
[1,]  1.0  0.5  -0.5
[2,]  0.5  1.0  0.5
[3,]  -0.5  0.5  1.0
> #方法二
> covdef(data)
     [,1] [,2] [,3]
[1,]  1.0  0.5  -0.5
[2,]  0.5  1.0  0.5
[3,]  -0.5  0.5  1.0
> cov(data)   
     [,1] [,2] [,3]
[1,]  1.0  0.5 -0.5
[2,]  0.5  1.0  0.5
[3,]  -0.5  0.5  1.0

2 多元正态随机数的产生

2.1 模拟产生多元正态随机数

library(MASS)
data <- matrix(c(1,2,3,2,1,3,3,1,2),ncol = 3,nrow  = 3)
n <- 10
#期望
means <- colMeans(data)    
#协方差矩阵
sigma <- covfun(data)        
#多元正太分布数据
mydata <- mvrnorm(n,means,sigma)              
mydata <- as.data.frame(mydata)
names(mydata)=c("x",'y','z')

结果显示
> mydata
       x         y           z
1   1.7841501 2.6839626  2.89981247
2   2.7186884 2.9585370  2.23984853
3  -1.3911512 0.4725601  3.86371129
4   1.2814312 1.1688948  1.88746361
5   2.4993147 1.7965445  1.29722984
6   2.3202226 3.8700157  3.54979296
7   1.8283135 1.2382355  1.40992200
8   0.4241805 1.7382137  3.31403322
9   3.8676690 1.8421394  -0.02552953
10  2.1399611 2.2350913  2.09513021

2.2 独立性的判断

由于 n n n 维正态分布随机向量 X = [ X 1 , X 2 , … , X n ] T \mathbf{X}=\left[X_{1}, X_{2}, \ldots, X_{n}\right]^{T} X=[X1​,X2​,…,Xn​]T 相互独立的充要条件是它们两两不相关。判断多维正态随机变量随机数之间是否独立可转换为**判断其协方差矩阵是否为对角矩阵。**下面给出相应的和结果：

#判断数据间的独立性，及判断协方差矩阵是否为对角矩阵
indepence <- cov(mydata)

结果显示
> indepence
      x         y           z
x  2.0347762  0.7886578   -1.246118
y  0.7886578  0.9641959   0.175538
z  -1.2461183  0.1755380   1.421656

结果分析：从多元正态随机数的协方差矩阵可以看出，其协方差矩阵为非对角矩阵，因此该组数据为不相互独立。

3 相关性的判断

3.1 相关性基础知识

假设数据集 T = { x i } i = 1 m T=\left\{x_{i}\right\}_{i=1}^{m} T={xi​}i=1m​表示m个样本，每个样本表示为 X i = ( x i 1 , x i 2 , … , x i n ) T {{X}_{i}}={{\left( {{x}_{i1}},{{x}_{i2}},\ldots ,{{x}_{in}} \right)}^{T}} Xi​=(xi1​,xi2​,…,xin​)T。所有样本可以组成一个m×n的矩阵。则样本的第 j j j个分量与第 k k k分量的相关系数为：

r j k = s j k s j j s k k , j , k = 1 , 2 , ⋯   , n {{r}_{jk}}=\frac{{{s}_{jk}}}{\sqrt{{{s}_{jj}}}\sqrt{{{s}_{kk}}}},\quad j,k=1,2,\cdots ,n rjk​=sjj​ ​skk​ ​sjk​​,j,k=1,2,⋯,n

称

R = [ r 11 r 12 ⋯ r 1 p r 21 r 22 ⋯ r 2 p ⋮ ⋮ ⋮ r p 1 r p 2 ⋯ r p p ] R=\left[\begin{array}{cccc}r_{11} & r_{12} & \cdots & r_{1 p} \\ r_{21} & r_{22} & \cdots & r_{2 p} \\ \vdots & \vdots & & \vdots \\ r_{p 1} & r_{p 2} & \cdots & r_{p p}\end{array}\right] R=⎣⎢⎢⎢⎡​r11​r21​⋮rp1​​r12​r22​⋮rp2​​⋯⋯⋯​r1p​r2p​⋮rpp​​⎦⎥⎥⎥⎤​

为样本的相关矩阵。

对于相关性检验有Pearson相关性检验，Spearman秩检验和Kendall秩检验。在R软件中，用函数cor.test() 进行检验。

若检验正态随机向量两两不相关，**利用正态分布的相关性与独立性等价，对相关的随机变量的协方差矩阵通过变化转化为协方差矩阵为对角矩阵。**下面给出相应的定理：

定理 若随机向量X 服从 N ( μ , B ) N(μ,B) N(μ,B),则存在一个正交变换 U U U,使得 Y = U X Y=UX Y=UX 是一个相互独立的正态随机向量。

证明： B为实对称矩阵, 存在正交阵 U U U, 使
U B U r = D = [ d 1 d 2 ⋱ d n ] \mathbf{U B U}^{r}=\mathbf{D}=\left[\begin{array}{llll} d_{1} & & & \\ & d_{2} & & \\ & & \ddots & \\ & & & d_{n} \end{array}\right] UBUr=D=⎣⎢⎢⎡​d1​​d2​​⋱​dn​​⎦⎥⎥⎤​

其中， d i d_{i} di​ 是 B B B 的特征值。又因 B B B是正定阵(从而非奇异的)，则 B B B 有 n n n个线性无关特征向量。

从而可得 U U U是以特征向量为列构成的正交阵。

最后通过 Y = U X Y=UX Y=UX 可将具有相关性的随机变量转化为相互独立的随机变量。

3.2 代码实现与结果分析

下面用函数cor() 和 cor.test() 求矩阵mydata的相关系数矩阵和相关性检验。下面给出相应的代码和结果。

#相关系数矩阵
> cor(mydata)
       x         y            z
x  1.0000000   0.5630509    -0.7326620
y  0.5630509   1.0000000    0.1499309
z  -0.7326620  0.1499309    1.0000000

cor.test(~x+y, mydata)
	Pearson's product-moment correlation

data:  x and y
t = 1.927, df = 8, p-value = 0.09013
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1031396  0.8805219
sample estimates:
      cor 
0.5630509 

> cor.test(~x+z, mydata)

	Pearson's product-moment correlation

data:  x and z
t = -3.0448, df = 8, p-value = 0.01595
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9322420 -0.1912684
sample estimates:
      cor 
-0.732662 

> cor.test(~y+z, mydata)

	Pearson's product-moment correlation

data:  y and z
t = 0.42892, df = 8, p-value = 0.6793
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.5296992  0.7123144
sample estimates:
      cor 
0.1499309

结果分析： 在 x x x和 y y y的相关系数检验中，其P值为0.09013>0.05,接受原假设，认为变量 x x x和 y y y不相关。在 x x x和 z z z的相关系数检验中，其P值为0.01595<0.05，拒绝原假设，认为变量 x x x和 z z z的具有相关性。在 y y y和 z z z的的相关系数检验中，其P值为0.6793>0.05，接受原假设，认为变量y和z不相关。

下面利用正态分布的相关性与独立性等价，对相关的随机变量的协方差矩阵通过变化转化为协方差矩阵为对角矩阵。以下是计算步骤：

Step1: 计算多元正态分布的协方差矩阵 B B B;

Step2: 计算协方差矩阵B的特征值和特征向量；

Step3: 以特征向量为列构成的正交阵 U U U;

Step4: 通过 Y = X U Y=XU Y=XU 将具有相关性的随机变量转化为相互独立的随机变量;

Step5: 计算经过变换后的随机数Y的协方差矩阵，判断是否为对角矩阵，从而检验结果的准确性。

下面给出相应的代码和结果。

ev <- eigen(cov(mydata))
#特征值
te <- ev$values
#特征向量
U <- as.matrix(ev$vectors)
U
Y <- matrix(nrow = n,ncol = 3)
Y <- as.matrix(mydata)%*%U
cov(Y)

> U
           [,1]      [,2]       [,3]
[1,]  0.7968303 0.1781240  0.5773503
[2,]  0.2441552 0.7791373 -0.5773503
[3,] -0.5526751 0.6010133  0.5773503
> Y <- matrix(nrow = n,ncol = 3)
> Y <- as.matrix(mydata)%*%U
> cov(Y)
        [,1]          [,2]          [,3]
[1,]  3.140725e+00 -1.304423e-16 -1.003642e-15
[2,] -1.304423e-16  1.279903e+00  1.216779e-15
[3,] -1.003642e-15  1.216779e-15  6.933583e-16

最后，从矩阵 Y Y Y的协方差矩阵可看出，其协方差矩阵为对角矩阵，从而把具有相关性的正态分布转化为了不相关的随机数。