第一种标准化转换公式:x*=D-1(x-µ),求出样本x的期望和其协方差矩阵的对角矩阵的逆即可。
第二种标准化转换公式:x*=∑-1/2(x-µ),其中∑-1/2=TΛT′,T为x的协方差矩阵的特征向量矩阵,Λ为x的协方差矩阵的特征值构成的对角矩阵。
第一种转换能够消除各变量单位的或方差差异的影响,但不能消除变量之间的相关性的影响。
第二种转换则可以做到消除变量之间的相关性的影响。
下面是R语言代码:
library(MASS) Sigma=matrix(c(2,2.2,2.2,7),2)#建立一个2行的协方差矩阵 mu=c(1,1)#期望 set.seed(8)#设置一个随机种子 x=mvrnorm(n=100,mu,Sigma)#随机生成100个随机样本服从期望为mu,协方差为Sigma的二元正态分布 s=cov(x)#计算X的协方差矩阵 plot(x,ylim=c(-10,10),xlim=c(-10,10),pch=20,ylab="y",xlab="x")#对随机样本x的点画图,横纵坐标数值范围为-10到10,画图符号选20号符号
#diag函数可以用来提取矩阵的对角元素,并将它保存成向量的格式,然后继续用diag函数来生成给定对角元素的对角阵 D=diag(sqrt(diag(s)))#计算对角矩阵D meanx=colMeans(x)#按列来求X的均值 #R语言是默认byrow=False,即数据按列输入,byrow=True则是按行输入,rep函数用于重复输出 mean=matrix(rep(meanx,100),nrow=100,byrow=T)#将meanx按行输入100次生成一个100*2的矩阵 x1=t(solve(D)%*%t(x-mean))#标准化转换,D是2*2矩阵,X-mean是100*2矩阵,所以要转置X-mean进行运算,最后再转置回来,%*%是左乘的意思 plot(x1,ylim=c(-10,10),xlim=c(-10,10),pch=20,ylab="y",xlab="x") cov(x1)#x1的协方差阵
完成第一种标准转换后,确实消除了消除x和y单位的或方差差异的影响,但显然没有消除x对y或y对x的相关性的影响。图中直观看到随着x增加y呈上升趋势,协方差阵cov(x,y),cov(y,x)不为0。下面进行第二种标准转换。
s.eigen=eigen(s)#计算s的特征值和特征向量 s.eigen$vectors#输出s的特征向量 a=s.eigen$values#s的特征值赋值给a A=diag(a)#计算对角矩阵 T=s.eigen$vectors#S的单位正交特征向量定义了一个新的坐标系 x2=t(t(T)%*%t(x))#将随机样本X投影到新的坐标系T上 plot(x2,ylim=c(-10,10),xlim=c(-10,10),pch=20,ylab="y",xlab="x") round(cov(x2),4)#求x2协方差,保留四位小数
在x落在由特征向量定义的坐标轴上,可以看到变量间的相关性已经被消除了。相当于将原先的坐标轴旋转了一下。下面就继续消除单位的或方差差异的影响
s1=T%*%sqrt(A)%*%t(T)#求A的平方根矩阵 x3=t(solve(s1)%*%t(x-mean))#标准化转换 plot(x3,ylim=c(-10,10),xlim=c(-10,10),pch=20,ylab="y",xlab="x") round(cov(x3),4)#求x3协方差,保留四位小数
最后再输出所有的plot图,直观的看看变换过程。