R多元统计分析上机讲义全.docx
《R多元统计分析上机讲义全.docx》由会员分享，可在线阅读，更多相关《R多元统计分析上机讲义全.docx（62页珍藏版）》请在冰豆网上搜索。
R多元统计分析上机讲义全
应用多元统计分析
应用多元统计分析
AppliedMultivariateStatisticalAnalysis
第一章绪论
在实际问题中，很多随机现象涉及到的变量不是一个，而是经常是多个变量，并且这些变量间又存在一定的联系。
我们经常需要处理多个变量的观测数据，如果用一元统计方法，由于忽视了各个变量之间可能存在的相关性，一般说来，丢失信息太多，分析的结果不能客观全面反映数据所包含的容，因此，我们就需要用到多元统计的方法。
多元统计分析（MultivariateStatisticalAnalysis）也称多变量统计分析、多因素统计分析或多元分析，是研究客观事物中多变量（多因素或多指标）之间的相互关系和多样品对象之间差异以及以多个变量为代表的多元随机变量之间的依赖和差异的现代统计分析理论和方法。
多元统计分析是解决实际问题的有效的数据处理方法。
随着电子计算机使用的日益普及，多元统计统计方法已广泛地应用于自然科学、社会科学的各个方面。
第二章矩阵
矩阵即是二维的数组，它非常的重要，以至于需要单独讨论。
由于矩阵应用非常广泛，因此对它定义了一些特殊的应用和操作，R包括许多只对矩阵操作的操作符和函数。
2.1矩阵的建立
在R中最为常用的是用命令matrix（）建立矩阵，而对角矩阵常用函数diag（）建立。
例如
>X<-matrix（1,nr=2,nc=2）
>X
[,1][,2]
[1,]11
[2,]11
>X<-diag（3）#生成单位阵
>X
[,1][,2][,3]
[1,]100
[2,]010
[3,]001
>diag（2.5,nr=3,nc=5）
[,1][,2][,3][,4][,5]
[1,]2.50.00.000
[2,]0.02.50.000
[3,]0.00.02.500
>X<-matrix（1:
4,2）#等价于X<-matrix（1:
4,2,2）
>X
[,1][,2]
[1,]13
[2,]24
>rownames（X）<-c（"a","b"）
>colnames（X）<-c（"c","d"）
>X
cd
a13
b24
>dim（X）
[1]22
>dimnames（X）
[[1]]
[1]"a""b"
[[2]]
[1]"c""d"
注意：
循环准则仍然适用于matrix（），但要求数据项的个数等于矩阵的列数的倍数,否则会出现警告。
矩阵的维数使用c（）会得到不同的结果（除非是方阵）,因此需要小心。
数据项填充矩阵的方向可通过参数byrow来指定,其缺省是按列填充的（byrow=FALSE），byrow=TRUE表示按行填充数据。
再看几个例子：
>X<-matrix（1:
4,2,4）#按列填充
>X
[,1][,2][,3][,4]
[1,]1313
[2,]2424
>X<-matrix（1:
4,2,3）
Warningmessage:
Inmatrix（1:
4,2,3）:
数据长度[4]不是矩阵列数[3]的整倍数
>X<-matrix（1:
4,c（2,3））#不经常使用
>X
[,1][,2]
[1,]13
[2,]24
>X<-matrix（1:
4,2,4,byrow=TRUE）#按行填充
>X
[,1][,2][,3][,4]
[1,]1234
[2,]1234
因为矩阵是数组的特例,R中数组由函数array（）建立,因此矩阵也可以用函数array（）来建立，其一般格式为：
>array（data,dim,dimnames）
其中data为一向量，其元素用于构建数组；dim为数组的维数向量（为数值型向量）；dimnames为由各维的名称构成的向量（为字符型向量），缺省为空。
看几个例子：
>A<-array（1:
6,c（2,3））
>A
[,1][,2][,3]
[1,]135
[2,]246
>A<-array（1:
4,c（2,3））
>A
[,1][,2][,3]
[1,]131
[2,]242
>A<-array（1:
8,c（2,3））
>A
[,1][,2][,3]
[1,]135
[2,]246
2.2矩阵的下标（index）与子集（元素）的提取
矩阵的下标可以使用正整数、负整数和逻辑表达式，从而实现子集的提取或修改。
考查矩阵
>x<-matrix（1:
6,2,3）
>x
[,1][,2][,3]
[1,]135
[2,]246
•提取一个元素
>x[2,2]
[1]4
•提取若一个或若干个行或列
>x[2,2]
[1]4
>x[2,]
[1]246
>x[,2]
[1]34
>x[,2,drop=FALSE]
[,1]
[1,]3
[2,]4
>x[,c（2,3）,drop=FALSE]
[,1][,2]
[1,]35
[2,]46
•去掉某一个或若干个行与列
>x[-1,]
[1]246
>x[,-2]
[,1][,2]
[1,]15
[2,]26
•添加与替换元素
>x[,3]<-NA
>x
[,1][,2][,3]
[1,]13NA
[2,]24NA
>x[is.na（x）]<-1#缺失值用1代替
>x
[,1][,2][,3]
[1,]131
[2,]241
2.3矩阵四则运算
矩阵也可以进行四则运算（“+”、“-”、“*”、“/”，“^”），分别解释为矩阵对应元素的四则运算。
在实际应用中，比较有实际应用的是矩阵的相加，相减，相乘和矩阵的求逆。
矩阵的加减运算一般要求矩阵形状完全相同（dim属性完全相同），矩阵的相乘一般要求一矩阵的列维数与另一矩阵的行维数相同，而矩阵要求逆的话，一般要求它为一方阵。
2.3.1矩阵的加减运算
若A，B为两个形状相同的矩阵，两矩阵的和为C，R中表达式为：
C<-A+B
两矩阵的差为D，R中表达式为：
D<-A-B
矩阵也可以与数进行加减，A+5表示A中的每个元素加上5。
2.3.2矩阵的相乘
操作符%*%用于矩阵相乘。
若矩阵A的列数等于矩阵B的行数，矩阵A乘以矩阵B表示为：
A%*%B
注：
X*Y表示两个矩阵的逐元相乘，而不是X和Y的乘积。
2.3.3矩阵的求逆
若矩阵A为一方阵，矩阵的逆可以用下面的命令计算：
solve（A）。
操作符solve（）可以用来求解线性方程组：
Ax=b，解为solve（A,b）
在数学上，用直接求逆的办法解x<-solve（A）%*%b相比solve（A,b）不仅低效而且还有一种潜在的不稳定性。
2.4矩阵的其他一些代数运算
2.4.1求转置矩阵
转置函数为t（），矩阵X的转置为t（X）。
2.4.2提取对角元素
提取对角元的函数为diag（）。
例如：
>X<-matrix（1:
4,2,2）
>diag（X）
[1]14
事实上，diag（）的作用依赖于自变量，diag（vector）返回以自变量（向量）为主对角元素的对角矩阵；diag（matrix）返回由矩阵的主对角元素所组成的向量；diag（k）（k为标量）返回k阶单位阵。
2.4.3矩阵的合并与拉直
函数cbind（）把几个矩阵横向拼成一个大矩阵，这些矩阵行数应该相同；函数rbind（）把几个矩阵列向拼成一个大矩阵，这些矩阵列数应该相同。
（如果参与合并的矩阵比其它矩阵行数少或列数少，则循环不足后合并。
）例如：
>m1<-matrix（1,nr=2,nc=2）
>m1
[,1][,2]
[1,]11
[2,]11
>m2<-matrix（2,nr=2,nc=2）
>m2
[,1][,2]
[1,]22
[2,]22
>rbind（m1,m2）
[,1][,2]
[1,]11
[2,]11
[3,]22
[4,]22
>cbind（m1,m2）
[,1][,2][,3][,4]
[1,]1122
[2,]1122
2.4.4方阵的行列式
求方阵的行列式使用det（）：
X<-matrix（1:
4,2）
>X
[,1][,2]
[1,]13
[2,]24
>det（X）
[1]-2
2.4.5矩阵的特征根和特征向量
函数eigen（）用来计算矩阵的特征值和特征向量。
这个函数的返回值是一个含有values和vectors两个分量的列表。
命令
A<-eigen（X）
>A
$values
[1]5.3722813-0.3722813
$vectors
[,1][,2]
[1,]-0.5657675-0.9093767
[2,]-0.82456480.4159736
2.4.6Matrixfacilites
Inthefollowingexamples,AandBarematricesandxandbareavectors.
OperatororFunction
Description
A*B
Element-wisemultiplication
A%*%B
Matrixmultiplication
A%o%B
Outerproduct.AB'
crossprod（A,B）
crossprod（A）
A'BandA'Arespectively.
t（A）
Transpose
diag（x）
Createsdiagonalmatrixwithelementsofxintheprincipaldiagonal
diag（A）
Returnsavectorcontainingtheelementsoftheprincipaldiagonal
diag（k）
Ifkisascalar,thiscreatesakxkidentitymatrix.Gofigure.
solve（A,b）
Returnsvectorxintheequationb=Ax（i.e.,A-1b）
solve（A）
InverseofAwhereAisasquarematrix.
ginv（A）
Moore-PenroseGeneralizedInverseofA.
ginv（A）requiresloadingtheMASSpackage.
y<-eigen（A）
y$valaretheeigenvaluesofA
y$vecaretheeigenvectorsofA
y<-svd（A）
SinglevaluedecompositionofA.
y$d=vectorcontainingthesingularvaluesofA
y$u=matrixwithcolumnscontaintheleftsingularvectorsofA
y$v=matrixwithcolumnscontaintherightsingularvectorsofA
R<-chol（A）
CholeskifactorizationofA.Returnstheuppertriangularfactor,suchthatR'R=A.
y<-qr（A）
QRdecompositionofA.
y$qrhasanuppertrianglethatcontainsthedecompositionandalowertrianglethatcontainsinformationontheQdecomposition.
y$rankistherankofA.
y$qrauxavectorwhichcontainsadditionalinformationonQ.
y$pivotcontainsinformationonthepivotingstrategyused.
cbind（A,B,...）
Combinematrices（vectors）horizontally.Returnsamatrix.
rbind（A,B,...）
Combinematrices（vectors）vertically.Returnsamatrix.
rowMeans（A）
Returnsvectorofrowmeans.
rowSums（A）
Returnsvectorofrowsums.
colMeans（A）
Returnsvectorofcolumnmeans.
colSums（A）
Returnsvectorofcoumnsums.
2.4.7.其它函数
交叉乘积（crossproduct）,函数为crossprod（），crossprod（X,Y）表示一般的积X′Y，即X的每一列与Y的每一列的积组成的矩阵；QR分解,函数为qr（）,矩阵X的QR分解为X=QR，Q为正交阵，R为上三角阵；等等。
2.5矩阵的统计运算
函数cov（）和cor（）分别用于计算矩阵的协方差阵和相关系数阵。
矩阵的排列是有方向性的，在R中规定矩阵是按列排的，若没有特别说明，函数max（），min（），median（），var（），sd（），sum（），cumsum（），cumprod（），cummax（），cummin（）的使用对于矩阵也是按列计算的，但也可以通过选项MARGIN来改变。
下面我们要用到对一个对象施加某种运算的函数apply（），其格式为
>apply（X,MARGIN,FUN）
其中X为参与运算的矩阵,FUN为上面的一个函数或“+”、“-”、“*”、“＼”（必须放在引号中），MARGIN=1表示按列计算，MARGIN=2表示按行计算。
我们还用到sweep（）函数，命令
>sweep（X,MARGIN,STATS,FUN）
表示从矩阵X中按MATGIN计算STATS，并从X中除去（sweepout）。
2.5.1求均值
>m<-matrix（rnorm（n=12）,nrow=3）
>apply（m,MARGIN=1,FUN=mean）#求各行的均值
[1]-0.37738650.38641380.2052353
>apply（m,MARGIN=2,FUN=mean）#求各列的均值
[1]0.33862020.7320669-0.4624578-0.3225460
2.5.2标准化
>scale（m,center=T,scale=T）
2.5.3减去中位数
>row.med<-apply（m,MARGIN=1,FUN=median）
>sweep（m,MARGIN=1,STATS=row.med,FUN=”-”）
第三章多元正态分布及参数的估计
3.1绘制二元正态密度函数及其相应等高线图
书上例2.2.2，
时的二元正态密度函数及其等高线图：
x<-seq（-3,3,by=0.1）
y<-x
f<-function（x,y,a=1,b=1,r=0）{
a1=sqrt（a）
b1=sqrt（b）
d=1-r*r
d1=sqrt（d）*a1*b1
z=1/（2*pi*d1）*exp（（-x*x/a-y*y/b+2*r*x*y/（a1*b1））/（2*d））
}
z<-outer（x,y,f）#外积函数
persp（x,y,z,xlim=range（x）,ylim=range（y）,zlim=range（z,na.rm=TRUE）,
theta=30,nticks=5,ticktype="detailed",sub="σ1=σ2=1，ρ=0时的二元正态密度函数"）
#密度函数图
contour（x,y,z）#等高线图
image（x,y,z）#等高线图，实际数据大小用不同色彩表示
所得图形为：
相应等高线图
Outer（x,y,f）是一个一般性的外积函数，调用函数f，把x的任一个元素与y的任意一个元素搭配起来作为f的自变量计算得到新的元素值，当函数缺省时表示乘积情况。
对参数进行修改，可以绘制任一二元正态密度函数及其相应的等高线图。
3.2多元正态分布的参数估计
3.2.1多元正态总体的相关量
设观测数据阵为
•样本均值向量
设
，
=1，2，
，
，则样本均值向量Xn：
由2.5.1可得：
>Xn<-apply（x,MARGIN=2,mean）
或者
>ln<-rep（1,n）
>Xn<-（ln%*%x）/n
Xn即为所求样本均值向量。
•样本离差阵（交叉乘积阵）
样本离差阵A：
。
>A<-crossprod（x）-2*Xn%*%t（Xn）
或者
>m<-diag（1,n）-matrix（1,n,n）/n
>A<-t（x）%*%m%*%x
A即为所求样本离差阵。
•样本协方差阵
R中求样本协方差阵的函数为cov（）。
样本数据阵X的协方差矩阵S即为：
>S<-cov（X）
•样本相关阵
R中求样本协方差阵的函数为cor（）。
样本数据阵X的协方差矩阵R即为:
>R<-cor（X）
3.2.2极大似然估计
极大似然估计法是建立在极大似然原理基础上的一种统计方法。
设总体X，其概率密度函数（连续情况）或分布律（离散情况）为
，其中
是未知参数（或未知参数向量）。
设X1，X2，…，Xn为取自总体X的样本，则似然函数
为：
…,
）=
求使似然函数达到最大的参数
的值，即极大似然估计值。
在单参数场合，在R中可以使用函数optimize（）求极大似然估计值。
optimize（）的调用格式如下：
optimize（f=,interval=,lower=min（interval）,
upper=max（interval）,maximum=TRUE,
tol=.Machine$double.eps^0.25,…）
说明：
f是似然函数，interval是参数
的取值围，lower是
的下界，upper是
的上界，maximum=TRUE是求极大值，否则（maximum=FALSE）表示求函数的极小值，tol是表示求值的精确度，…是对f的附加说明。
在多参数场合，在R中用函数optim（）或者nlm（）来求似然函数的极大值，并求相应的极大值点。
optim（）的调用格式如下：
optim（par,fn,gr=NULL,
method=c（"Nelder-Mead","BFGS","CG","L-BFGS-B","SANN"）,
lower=-Inf,upper=Inf,
control=list（）,hessian=FALSE,…）
nlm（）的定义如下：
nlm（f,p,hessian=FALSE,typsize=rep（1,length（p））,fscale=1,
print.level=0,ndigit=12,gradtol=1e-6,
stepmax=max（1000*sqrt（sum（（p/typsize）^2））,1000）,
steptol=1e-6,iterlim=100,check.analyticals=TRUE,…）
三者主要区别是：
函数nlm（）仅使用牛顿-拉夫逊算法求函数的最小值点；函数optim（）提供method选项给出的5种方法中的一种进行优化；上面二个可用于多维函数的极值问题,，而函数optimize（）仅适用于一维函数，但可以用于最大与最小值点。
（具体选项见帮助。
）
第四章多元正态总体参数的假设检验
在一元统计中，用于检验一元正态总体参数
，
的抽样分布有
分布，
分布、F分布风，它们都是来自总体
的随机样本导出的检验统计量。
推广到多元正态总体后，也有相应于以上三个常用分布的统计量：
威沙特（Wishart）统计量，霍特林（Hotelling）
统计量，威尔克斯（Wilks）
统计量，这些统计量是多元统计分析所涉及的假设检验问题的基础。
4.1几个重要统计量的分布
对于多元正态总体来说,存在几个重要的统计量:
威沙特（Wishart）统计量，霍特林（Hotelling）
统计量，威尔克斯（Wilks）
统计量等,讨论这些统计量的分布是多元统计分析所涉及的假设检验问题的基础。
4.2单总体均值向量的检验及置信域
4.2.1均值向量的检验
书上例3.2.1，R程序如下
>x<-matrix（c（3.7,48.5,9.3,5.7,65.1,8.0,3.8,47.2,10.9,3.2,53.2,12.0,
3.1,55.5,9.7,4.6,36.1,7.9,2.4,24.8,14.0,7.2,33.1,7.6,6.7,47.4,8.5,
5.4,54.1,11.3,3.9,36.9,12.7,4.5,58.8,12.3,3.5,27.8,9.8,4.5,40.2,8.4,
1.5,13.5,10.1,8.5,56.4,7.1,4.5,71.6,8.2,6.5,52.8,10.9,4.1,44.1,11.2,
5.5,40.9,9.4）,20,3,byrow=TRUE）
>n<-20
>p<-3
>u0<-c（4,50,10）#所给总体均值
>ln<-rep（1,20）
>x0<-（ln%*%x）/n#样本均值
>xm<-x0-u0
>mm<-diag（1,20）-matrix（1,20,20）/n
>a<-t（x）%*%mm%*%x#样本离差阵
>ai=solve（a）
>dd=xm%*%ai%*%t（xm）
>d2=（n-1）*dd
>t2=n*d2;
>f<-（n-p）*t2/（（n-1）*p）#检验统计量
>f
[,1]
[1,]2.904546
>fa<-qf（0.95,p,n-p）#自由度为（p,n-p）的F分布的0
升级会员 