R-多元统计分析上机讲义

应用多元统计分析R实验上机讲义
应用多元统计分析 (4)
Applied Multivariate Statistical Analysis (4)
第一章绪论 (4)
第二章矩阵 (4)
2.1矩阵的建立 (4)
2.2矩阵的下标(index)与子集(元素)的提取 (6)
2.3 矩阵四则运算 (7)
2.3.1 矩阵的加减运算 (7)
2.3.2 矩阵的相乘 (8)
2.3.3 矩阵的求逆 (8)
2.4矩阵的其他一些代数运算 (8)
2.4.1 求转置矩阵 (8)
2.4.2 提取对角元素 (8)
2.4.3矩阵的合并与拉直 (8)
2.4.4方阵的行列式 (9)
2.4.5 矩阵的特征根和特征向量 (9)
2.4.6 其它函数 (9)
2.5 矩阵的统计运算 (11)
2.5.1 求均值 (11)
2.5.2 标准化 (11)
2.5.3 减去中位数 (11)
第三章多元正态分布及参数的估计 (12)
3.1 绘制二元正态密度函数及其相应等高线图 (12)
3.2 多元正态分布的参数估计 (14)
3.2.1 多元正态总体的相关量 (14)
3.2.2 极大似然估计 (14)
第四章多元正态总体参数的假设检验 (15)
4.1 几个重要统计量的分布 (15)
4.2 单总体均值向量的检验及置信域 (16)
4.2.1均值向量的检验 (16)
4.2.2样本协方差阵的特征值和特征向量 (17)
4.3多总体均值向量的检验 (17)
4.3.1 两正态总体均值向量的检验 (17)
4.3.2 多个正态总体均值向量的检验-多元方差分析 (19)
4.4协方差阵的检验 (20)
4.4.2 多总体协方差阵的检验 (20)
4.5独立性检验 (20)
4.6正态性检验 (21)
第五章判别分析 (22)
5.1距离判别 (22)
5.1.1 马氏距离 (22)
5.1.2 两总体的距离判别 (22)
5.1.3 多个总体的距离判别 (26)
5.2贝叶斯判别法及广义平方距离判别法 (26)
5.2.1 先验概率(先知知识) (26)
5.2.2 广义平方距离 (26)
5.2.3 后验概率（条件概率） (27)
5.2.4 贝叶斯判别准则 (27)
5.3费希尔（Fisher）判别 (29)
第六章聚类分析 (30)
6.2距离和相似系数 (30)
6.2.1距离 (31)
6.2.2数据中心化与标准化变换 (31)
6.2.3相似系数 (31)
6.3 系统聚类法 (31)
6.4类个数的确定 (34)
6.5动态聚类法 (36)
6.7变量聚类方法 (36)
第七章主成分分析 (37)
7.2 样本的主成分 (38)
7.3 主成分分析的应用 (39)
第八章因子分析 (42)
8.3 参数估计方法 (42)
8.4 方差最大的正交旋转 (45)
8.5 因子得分 (45)
第九章对应分析方法 (46)
第十章典型相关分析 (48)
应用多元统计分析
Applied Multivariate Statistical Analysis
第一章绪论
在实际问题中，很多随机现象涉及到的变量不是一个，而是经常是多个变量，并且这些变量间又存在一定的联系。

我们经常需要处理多个变量的观测数据，如果用一元统计方法，由于忽视了各个变量之间可能存在的相关性，一般说来，丢失信息太多，分析的结果不能客观全面反映数据所包含的内容，因此，我们就需要用到多元统计的方法。

多元统计分析(Multivariate Statistical Analysis)也称多变量统计分析、多因素统计分析或多元分析，是研究客观事物中多变量(多因素或多指标)之间的相互关系和多样品对象之间差异以及以多个变量为代表的多元随机变量之间的依赖和差异的现代统计分析理论和方法。

多元统计分析是解决实际问题的有效的数据处理方法。

随着电子计算机使用的日益普及，多元统计统计方法已广泛地应用于自然科学、社会科学的各个方面。

第二章矩阵
矩阵即是二维的数组，它非常的重要，以至于需要单独讨论。

由于矩阵应用非常广泛，因此对它定义了一些特殊的应用和操作，R 包括许多只对矩阵操作的操作符和函数。

2.1矩阵的建立
在R中最为常用的是用命令matrix( )建立矩阵，而对角矩阵常用函数diag( )建立。

例如
> X<-matrix(1,nr=2,nc=2)
> X
[,1] [,2]
[1,] 1 1
[2,] 1 1
> X <- diag(3) # 生成单位阵
> X
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> diag(2.5, nr = 3, nc = 5)
[,1] [,2] [,3] [,4] [,5]
[1,] 2.5 0.0 0.0 0 0
[2,] 0.0 2.5 0.0 0 0
[3,] 0.0 0.0 2.5 0 0
> X <- matrix(1:4, 2) # 等价于X <- matrix(1:4, 2, 2)
> X
[,1] [,2]
[1,] 1 3
[2,] 2 4
> rownames(X) <- c("a", "b")
> colnames(X) <- c("c", "d")
> X
c d
a 1 3
b 2 4
> dim(X)
[1] 2 2
> 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,] 1 3 1 3
[2,] 2 4 2 4
> X <- matrix(1:4, 2, 3)
Warning message:
In matrix(1:4, 2, 3) : 数据长度[4]不是矩阵列数[3]的整倍数
> X <- matrix(1:4, c(2, 3)) # 不经常使用
> X
[,1] [,2]
[1,] 1 3
[2,] 2 4
> X <- matrix(1:4, 2, 4, byrow=TRUE) # 按行填充
> X
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 1 2 3 4
因为矩阵是数组的特例,R中数组由函数array( )建立, 因此矩阵也可以用函数array( )来建立，其一般格式为：
> array(data, dim, dimnames)
其中data为一向量，其元素用于构建数组；dim为数组的维数向量（为数值型向量）；dimnames为由各维的名称构成的向量（为字符型向量)，缺省为空。

看几个例子：
> A<-array(1:6,c(2,3))
> A
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> A<-array(1:4,c(2,3))
> A
[,1] [,2] [,3]
[1,] 1 3 1
[2,] 2 4 2
> A<-array(1:8,c(2,3))
> A
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
2.2矩阵的下标(index)与子集(元素)的提取
矩阵的下标可以使用正整数、负整数和逻辑表达式，从而实现子集的提取或修改。

考查矩阵
> x <- matrix(1:6, 2, 3)
> x
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
•提取一个元素
> x[2,2]
[1] 4
•提取若一个或若干个行或列
> x[2,2]
[1] 4
> x[2,]
[1] 2 4 6
> x[,2]
[1] 3 4
> x[,2,drop=FALSE]
[,1]
[1,] 3
[2,] 4
> x[,c(2,3),drop=FALSE]
[,1] [,2]
[1,] 3 5
[2,] 4 6
•去掉某一个或若干个行与列
> x[-1,]
[1] 2 4 6
> x[,-2]
[,1] [,2]
[1,] 1 5
[2,] 2 6
•添加与替换元素
> x[,3] <- NA
> x
[,1] [,2] [,3]
[1,] 1 3 NA
[2,] 2 4 NA
> x[is.na(x)] <- 1 # 缺失值用1代替
> x
[,1] [,2] [,3]
[1,] 1 3 1
[2,] 2 4 1
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] 1 4
事实上，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,] 1 1
[2,] 1 1
> m2 <- matrix(2, nr = 2, nc = 2)
> m2
[,1] [,2]
[1,] 2 2
[2,] 2 2
> rbind(m1, m2)
[,1] [,2]
[1,] 1 1
[2,] 1 1
[3,] 2 2
[4,] 2 2
> cbind(m1, m2)
[,1] [,2] [,3] [,4]
[1,] 1 1 2 2
[2,] 1 1 2 2
2.4.4方阵的行列式
求方阵的行列式使用det( )：
X<-matrix(1:4, 2)
> X
[,1] [,2]
[1,] 1 3
[2,] 2 4
> 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.8245648 0.4159736
2.4.6 Matrix facilites
In the following examples, A and B are matrices and x and b are a vectors.
Operator or Function
Description
A * B
Element-wise multiplication
A %*% B
Matrix multiplication
A %o% B
Outer product. AB'
crossprod(A,B)
crossprod(A)
A'B and A'A respectively.
t(A)
Transpose
diag(x)
Creates diagonal matrix with elements of x in the principal diagonal
diag(A)
Returns a vector containing the elements of the principal diagonal
diag(k)
If k is a scalar, this creates a k x k identity matrix. Go figure.
solve(A, b)
Returns vector x in the equation b = Ax (i.e., A-1b)
solve(A)
Inverse of A where A is a square matrix.
ginv(A)
Moore-Penrose Generalized Inverse of A.
ginv(A) requires loading the MASS package.
y<-eigen(A)
y$val are the eigenvalues of A
y$vec are the eigenvectors of A
y<-svd(A)
Single value decomposition of A.
y$d = vector containing the singular values of A
y$u = matrix with columns contain the left singular vectors of A
y$v = matrix with columns contain the right singular vectors of A
R <- chol(A)
Choleski factorization of A. Returns the upper triangular factor, such that R'R = A.
y <- qr(A)
QR decomposition of A.
y$qr has an upper triangle that contains the decomposition and a lower triangle that contains information on the Q decomposition.
y$rank is the rank of A.
y$qraux a vector which contains additional information on Q.
y$pivot contains information on the pivoting strategy used.
cbind(A,B,...)
Combine matrices(vectors) horizontally. Returns a matrix.
rbind(A,B,...)
Combine matrices(vectors) vertically. Returns a matrix.
rowMeans(A)
Returns vector of row means.
rowSums(A)
Returns vector of row sums.
colMeans(A)
Returns vector of column means.
colSums(A)
Returns vector of coumn sums.
2.4.7.其它函数
交叉乘积(cross product), 函数为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中除去(sweep out)。

2.5.1 求均值
> m<-matrix(rnorm(n=12),nrow=3)
> apply(m, MARGIN=1, FUN=mean) # 求各行的均值
[1] -0.3773865 0.3864138 0.2052353
> apply(m, MARGIN=2, FUN=mean) # 求各列的均值
[1] 0.3386202 0.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，…，X n为取自总体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,3
6.9,12.7,4.5,58.8,12.3,3.5,2
7.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,5
2.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.95分位数
> fa
[1] 3.196777
> b<-1-pf(f,p,n-p) # 尾概率值
> b
[,1]
[1,] 0.06492834
> beta<-pf(fa,p,n-p,t2) # 犯第二类错误的概率（假设总体均值）
> beta
[1] 0.3616381
取检验水平为0.05，由尾概率值p=0.064928340.05=，可得相容；同样由
F=2.904546 3.196777=Fa，也可得相容。

在这种情况下，可能犯第二类错误，概率为
=0.3616（假定总体均值）。

4.2.2样本协方差阵的特征值和特征向量
书上例3.2.2，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,3
6.9,12.7,4.5,58.8,12.3,3.5,2
7.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,5
2.8,10.9,4.1,44.1,11.2,
5.5,40.9,9.4),20,3,byrow=TRUE)
> s<-cov(x)
> s
[,1] [,2] [,3]
[1,] 2.879368 10.0100 -1.809053
[2,] 10.010000 199.7884 -5.640000
[3,] -1.809053 -5.6400 3.627658
> a<-eigen(s)
> a
$values
[1] 200.462464 4.531591 1.301392
$vectors
[,1] [,2] [,3]
[1,] -0.05084144 -0.57370364 0.81748351
[2,] -0.99828352 0.05302042 -0.02487655
[3,] 0.02907156 0.81734508 0.57541452
4.3多总体均值向量的检验
4.3.1 两正态总体均值向量的检验
书上例3.3.1，R程序为：
> n<-10
> m<-10
> p<-4
> x<-matrix(c(65,75,60,75,70,55,60,65,60,55,
+ 35,50,45,40,30,40,45,40,50,55,25,20,35,40,30,
+ 35,30,25,30,35,60,55,65,70,50,65,60,60,70,75),10)
> ln<-rep(1,n)
> x0<-(ln%*%x)/n
> mx<-diag(1,n)-matrix(1,n,n)/n
> a1<-t(x)%*%mx%*%x
> y<-matrix(c(55,50,45,50,55,60,65,50,40,45,
+ 55,60,45,50,50,40,55,60,45,50,40,45,35,50,
+ 30,45,45,35,30,45,65,70,75,70,75,60,75,80,65,70),10) > y0<-(ln%*%y)/n
> my<-diag(1,n)-matrix(1,n,n)/n
> a2<-t(y)%*%my%*%y
> a<-a1+a2
> xy<-x0-y0
> ai<-solve(a)
> dd<-xy%*%ai%*%t(xy)
> d2<-(m+n-2)*dd
> t2<-n*m*d2/(n+m)
> f<-(n+m-1-p)*t2/((n+m-2)*p)
> pp<-1-pf(f,p,m+n-p-1)
> x0
[,1] [,2] [,3] [,4]
[1,] 64 43 30.5 63
> y0
[,1] [,2] [,3] [,4]
[1,] 51.5 51 40 70.5
> a1
[,1] [,2] [,3] [,4]
[1,] 490 -170 -120.0 -245
[2,] -170 510 10.0 310
[3,] -120 10 322.5 260
[4,] -245 310 260.0 510
> a2
[,1] [,2] [,3] [,4]
[1,] 502.5 60 175 -7.5
[2,] 60.0 390 50 195.0
[3,] 175.0 50 450 -100.0
[4,] -7.5 195 -100 322.5
> d2
[,1]
[1,] 5.972499
> t2
[,1]
[1,] 29.86250
> f
[,1]
[1,] 6.221353
> pp
[,1]
[1,] 0.003705807
取检验水平为0.01，根据尾概率值p=0.0037058070.01=，可得应否定。

4.3.2 多个正态总体均值向量的检验-多元方差分析
书上例3.3.2，可利用类似例3.2.1或例3.3.1的程序进行计算得出结论。

下面我们用R自带的manova()函数进行分析。

程序如下：
x<-read.table("D:/data/d332.txt",header=T)
x<-as.matrix(x[,1:4])
rate<-factor(gl(3,20),labels=c("group1","group2","group3"))
fit <- manova(x ~ rate)
summary.aov(fit) # 对每一个变量进行单因素方差分析
summary(fit, test="Wilks") # 使用威尔克斯统计量
程序结果：
> summary.aov(fit)
Response x1 :
Df Sum Sq Mean Sq F value Pr(>F)
rate 2 39066 19533 8.878 0.0004401 ***
Residuals 57 125409 2200
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response x2 :
Df Sum Sq Mean Sq F value Pr(>F)
rate 2 4017 2009 2.8293 0.06738 .
Residuals 57 40467 710
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response x3 :
Df Sum Sq Mean Sq F value Pr(>F)
rate 2 13.43 6.72 0.1838 0.8326
Residuals 57 2082.50 36.54
Response x4 :
Df Sum Sq Mean Sq F value Pr(>F)
rate 2 17.20 8.60 0.4785 0.6222
Residuals 57 1024.40 17.97
> summary(fit, test="Wilks")
Df Wilks approx F num Df den Df Pr(>F)
rate 2 0.66212 3.09069 8 108 0.003538 **
Residuals 57
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
结果说明：
（1）取检验水平为0.01，则对四个指标逐项用一元方差分析方法进行检验，由p
值可得三个组指标间只有第一个指标有显著差异(=0.0004401)；
（2）取检验水平为0.01，利用威尔克斯统计量得到p=0.0035380.01，故拒绝
原假设，即认为三个组的指标之间有显著差异。

4.4协方差阵的检验
4.4.2 多总体协方差阵的检验
书上例3.4.1，R程序略（类似例3.2.1或例3.3.1）
4.5独立性检验
书中例3.5.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
> x0<-(ln%*%x)/n
> xm<-x0-u0
> mm<-diag(1,20)-matrix(1,20,20)/n
> a<-t(x)%*%mm%*%x
> a0<-det(a)
> a1<-a[1,1]*a[2,2]*a[3,3]
> v<-a0/a1
> b<-n-1.5-(p*p*p-3)/(3*p*p-3*3)
> df<-0.5*(p*(p+1)-2*3)
> kc<--b*log(v)
> p0<-1-pchisq(kc,df)
> kc
[1] 9.755514
> p0
[1] 0.02076288
取检验水平为0.05，根据尾概率值p=0.020762880.05=，可得应否定，由R软件所的结果与SAS软件所的结果一致。

4.6正态性检验
书中例3.6.1，R程序为：
> x<-matrix(c(100,99,96,99,96,75,97,68,76,62,67,34,100,97,100,96,78,97, + 89,88,84,39,78,37),12)
> n<-12
> p<-2
> ln<-rep(1,n)
> x0<-(ln%*%x)/n
> s<-cov(x)
> si<-solve(s)
> m<-0
> for(i in 1:n){
+ xx0<-x[i,]-x0
+ dd<-xx0%*%si%*%t(xx0)
+ print(c(i,dd))
+ if(dd<=1.386) m<-m+1
+ }
[1] 1.000000 0.883192
[1] 2.0000000 0.7787306
[1] 3.000000 0.696518
[1] 4.000000 0.789136
[1] 5.000000 2.188154
[1] 6.000000 2.384875
[1] 7.0000000 0.8767929
[1] 8.000000 2.033652
[1] 9.0000000 0.2691041
[1] 10.000000 5.046531
[1] 11.0000000 0.7891688
[1] 12.000000 5.264147
> m
[1] 7
> pp<-m/n
> pp
[1] 0.5833333
第五章判别分析
判别分析是用于判断样品所属类型的一种统计分析方法。

判别分析的目的是对已知归类的数据建立由数值指标构成的归类规则，然后把这样的规则应用到未知归类的样品去归类。

在生产、科研和日常生活中经常会遇到如何根据观测到的数据资料对所研究的对象进行判别归类的问题。

判别分析问题一般可以如下描述：设有k个维总体，其分布特征已知（如已知分布函数分别为，，或知道来自各个总体的训练样本）。

对给定的一个新样品X，判断它来自哪个总体。

通常我们先对预先得到的来自这k个总体的若干个样品(称为训练样品)进行检验和归类，来决定相应的判别归类问题是否有意义及误判可能性大小。

然后再对给定的一个或几个新的样品，进行判别归类，即决定它(们)自哪个总体。

解决这个问题可以有多种途径, 下面我们分别讨论几种常用的方法，如距离判别、贝叶斯判别、Fisher判别等。

R通用程序：
首先我们要用命令
>library(MASS)
MASS包里的lda( )针对线性判别分析。

加载MASS宏包，再用函数lda( )就可完成判别分析，其基本调用格式如下：
lda(formula, data, ... , subset, na.action)
说明：formula用法为groups~x1 + x2 +…，group表明总体来源，x1，x2，…表示分类指标；subset指明训练样本。

具体说明见R帮助。

5.1距离判别
5.1.1 马氏距离
马氏距离定义：样本X和总体
其中为总体均值向量，为总体协方差阵。

5.1.2 两总体的距离判别
判别准则：
其中=，为X到总体的距离。

关于两总体距离判别的R程序(参考薛毅教授的《统计建模与R软件》一书)：
discriminiant.distance <- function
(TrnX1, TrnX2, TstX = NULL, var.equal = FALSE){
if (is.null(TstX) == TRUE) TstX <- rbind(TrnX1,TrnX2)
if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
else if (is.matrix(TstX) != TRUE)
TstX <- as.matrix(TstX)
if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)
if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)
nx <- nrow(TstX)
blong <- matrix(rep(0, nx), nrow=1, byrow=TRUE,
dimnames=list("blong", 1:nx))
mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)
if (var.equal == TRUE || var.equal == T){
S <- var(rbind(TrnX1,TrnX2))
w <- mahalanobis(TstX, mu2, S) - mahalanobis(TstX, mu1, S)
}
else{
S1 < -var(TrnX1); S2 <- var(TrnX2)
w <- mahalanobis(TstX, mu2, S2) - mahalanobis(TstX, mu1, S1)
}
for (i in 1:nx){
if (w[i] > 0)
blong[i] <- 1
else
blong[i] <- 2
}
blong
}
在程序中，输入变量TrnX1、TrnX2表示训练样本X1，X2，其输入格式是数据框，或矩阵（样本按行输入），输入变量TstX是待测样本，其输入格式是数据框，或矩阵（样本按行输入），或向量（一个待测样本）。

如果不输入TstX（缺省值），则待测样本为两个训练样本之和，即训练样本的回代情况。

输入变量var.equal是逻辑变量，var.equal == TRUE表示两个总体的协方差相同；否则（缺省值）为不同。

在上述程序中，用到马氏距离函数mahalanobis()，该函数的使用格式为
mahalanobis (x, center, cov, inverted=FALSE, ...)
其中x是样本数据构成的向量或矩阵（p维），center为样本中心，cov为样本的协方差阵。

对于书中例5.1.1，调用discriminiant.distance( )进行判别（假设协方差阵相等）：
x<-data.frame(x1=c(13.85,22.31,28.82,15.29,28.79),
x2=c(2.79,4.67,4.63,3.54,4.90),
x3=c(7.8,12.31,16.18,7.50,16.12),
x4=c(49.6,47.8,62.15,43.20,58.10))
y<- data.frame(x1=c(2.18,3.85,11.40,3.66,12.10), x2=c(1.06,0.80,0.00,2.42,0.00),
x3=c(1.22,4.06,3.50,2.14,5.68),
x4=c(20.60,47.10,0.00,15.10,0.00))
testx<- rbind(
c(8.85, 3.38, 5.17, 26.10), c(28.60, 2.40, 1.20, 127.0),
c(20.70, 6.70, 7.60, 30.20), c(7.90, 2.40, 4.30, 33.20),
c(3.19, 3.20, 1.43, 9.90), c(12.40, 5.10, 4.43, 24.60),
c(16.80, 3.40, 2.31, 31.30), c(15.00, 2.70, 5.02, 64.00))
discriminiant.distance(x, y,var.equal=TRUE)
1 2 3 4 5 6 7 8 9 10
blong 1 1 1 1 1 2 2 2 2 2
discriminiant.distance(x, y,testx,var.equal=TRUE)
1 2 3 4 5 6 7 8
blong 2 1 1 2 2 1 1 1
由程序结果可得待判样品2，3，6，7，8属于含钾盐泉(A盆地)，其余三个属于不含钾盐泉(B盆地)。

利用R自带函数lda()，该例的R程序如下：
> w <- read.table(file="D:/data/disc511.txt",header=T)
> attach(w)
> names(w)
> library(MASS)
> z <- lda(group~x1+x2+x3+x4)
> z
> pred<-predict(z) $ class # predict( )是R内置函数，可以将lda( )的输出应用
于训练样品数据进行预测，从而进行对比。

> table(pred,group)
> newdata<-rbind(
+ c(8.85, 3.38, 5.17, 26.10), c(28.60, 2.40, 1.20, 127.0),
+ c(20.70, 6.70, 7.60, 30.20), c(7.90, 2.40, 4.30, 33.20),
+ c(3.19, 3.20, 1.43, 9.90), c(12.40, 5.10, 4.43, 24.60),
+ c(16.80, 3.40, 2.31, 31.30), c(15.00, 2.70, 5.02, 64.00))
> newdata<-data.frame(newdata)
> predict(z, newdata=newdata)
> detach(w)
R程序结果:
Call:
lda(group ~ x1 + x2 + x3 + x4)
Prior probabilities of groups:
A B
0.5 0.5
Group means:
x1 x2 x3 x4
A 21.812 4.106 11.982 52.17
B 6.638 0.856 3.320 16.56
Coefficients of linear discriminants: LD1
x1 -0.7794490
x2 -0.6888651
x3 1.4115135
x4 -0.1192217
group
pred A B
A 5 0
B 0 5
$class
[1] B A A B B A A A
Levels: A B
$posterior
A B
1 1.639701e-03 9.983603e-01
2 1.000000e+00 1.932625e-83
3 1.000000e+00 1.269619e-20
4 8.302424e-02 9.169758e-01
5 1.190922e-0
6 9.999988e-01
6 1.000000e+00 1.129611e-10
7 1.000000e+00 1.161894e-26
8 1.000000e+00 7.135903e-22
$x
LD1
1 1.0536512
2 -31.2985593
3 -7.5286829
4 0.3947245
5 2.2416596
6 -3.7639282
7 -9.8136273
8 -8.0017623
结果说明：
1) Group means: 包含了每组的平均向量；
2) Coefficients of linear discriminants: 线性判别系数；
3) 列联表表明将训练样品数据代入线性判别函数后的判别结果，两组都没有错判；
4) 由$class可以看出8个待判样品，待判样品2，3，6，7，8属于含钾盐泉(A盆地)，其余三个属于不含钾盐泉(B盆地)（与上一程序结果一致）；
5) $posterior给出了后验概率值（具体概念见5.2节）；
6)$x给出了线性判别函数的数值。

5.1.3 多个总体的距离判别
类似与两个总体的情况，多个总体的情况，按照距离最近的原则对X进行判别归类时，首先计算样品到各类的马氏(Mahalanobis)距离，然后进行比较，把待判样品判归距离最小的那个总体。

（自编关于多个总体距离判别的R函数可参考《统计建模与R软件一书》）。

5.2贝叶斯判别法及广义平方距离判别法
5.2.1 先验概率(先知知识)
设有k个总体，假设事先对所研究的问题有一定的认识，这种认识常用先验概率来描述，即已知这k个总体各自出现的概率（验前概率）为（显然，=1），这组验前概率称为先验概率。

5.2.2 广义平方距离
样品X到总体(=1,…,k)的广义平方距离为：
,
其中
是样品X到总体的马氏距离；
其中为第类的组内样本协方差阵。

5.2.3 后验概率（条件概率）
当样品X已知时，它属于的概率就称为后验概率，一般记为(或)。

5.2.4 贝叶斯判别准则
几个概念：
1.错判概率和错判损失；
2.关于先验概率的平均损失。

定义5.2.1：
设有k个总体：，相应的先验概率为（，
=1）。

如果有判别法，使得带来的平均损失达最小，则称判别
法符合贝叶斯判别准则，或称为贝叶斯判别的解。

出于例题需要，学习多总体的Bayes判别程序（两总体情况参考《统计建模与R软件一书》）：
distinguish.bayes <- function
(TrnX, TrnG, p = rep(1, length(levels(TrnG))),
TstX = NULL, var.equal = FALSE){
if ( is.factor(TrnG) == FALSE){
mx <- nrow(TrnX); mg <- nrow(TrnG)
TrnX <- rbind(TrnX, TrnG)
TrnG <- factor(rep(1:2, c(mx, mg)))
}
if (is.null(TstX) == TRUE) TstX <- TrnX
if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
else if (is.matrix(TstX) != TRUE) TstX <- as.matrix(TstX)
if (is.matrix(TrnX) != TRUE) TrnX <- as.matrix(TrnX)
nx <- nrow(TstX)
blong <- matrix(rep(0, nx), nrow=1, dimnames=list("blong", 1:nx))
g <- length(levels(TrnG))
mu <- matrix(0, nrow=g, ncol=ncol(TrnX))
for (i in 1:g)
mu[i,] <- colMeans(TrnX[TrnG==i,])
D <- matrix(0, nrow=g, ncol=nx)
if (var.equal == TRUE || var.equal == T){
for (i in 1:g){
d2 <- mahalanobis(TstX, mu[i,], var(TrnX))
D[i,] <- d2 - 2*log(p[i])
}
}
else{
for (i in 1:g){
S <- var(TrnX[TrnG==i,])
d2 <- mahalanobis(TstX, mu[i,], S)
D[i,] <- d2 - 2*log(p[i])-log(det(S))
}
}
for (j in 1:nx){
dmin <- Inf
for (i in 1:g)
if (D[i,j] < dmin){
dmin <- D[i,j]; blong[j] <- i
}
}
blong
}
在程序中，相同变量意义与距离判别程序中相同。

TrnG是因子变量，表示训练样本的分
类情况。

输入变量p是先验概率，缺省值均为1.
书上例5.2.2：
（1）调用函数distinguish.bayes ()进行判别（假设协方差阵不相等）：
X<-data.frame(x1=c(228,245,200,170,100,225,130,150,120,160,
185,170,165,135,100),
x2=c(134,134,167,150,167,125,100,117,133,100,115,125,142,108,117), x3=c(20,10,12,7,20,7,6,7,10,5,5,6,5,2,7),
x4=c(11,40,27,8,14,14,12,6,26,10,19,4,3,12,2))
v <- read.table(file="D:/data/d522.txt",header=T)
X<-v[,2:5]
G<-gl(3,5)
distinguish.bayes(X,G)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
blong 1 1 1 1 1 1 3 3 2 2 3 3 1 3 3
所得结果与书本不一致.
（2）利用R自带函数qda()进行二次判别分析，该例的R程序如下：
v <- read.table(file="D:/data/d522.txt",header=T)
attach(v)
names(v)
library(MASS)
qda.v<-qda(group~x1+x2+x3+x4) # qda()用于二次判别分析qda.v
pred.v<-predict(qda.v) $ class
table(pred.v,group)
R程序结果:
Call:
qda(group ~ x1 + x2 + x3 + x4)
Prior probabilities of groups:
1 2 3
0.3333333 0.3333333 0.3333333
Group means:
x1 x2 x3 x4
1 188.6 150.4 13.8 20.0
2 157.0 115.0 7.0 13.6
3 151.0 121.
4 5.0 8.0
group
pred.v 1 2 3
1 5 0 0
2 0 5 0
3 0 0 5
由结果可得通过二次判别分析，三组都没有错判。

5.3费希尔（Fisher）判别
书上例5.3.2,利用lda()函数进行费希尔判别的R程序如下：
> v <- read.table(file="D:/data/d522.txt",header=T)
> attach(v)
> names(v)
> library(MASS)
> v.lda<-qda(group~x1+x2+x3+x4)
> v.lda
> v.pred<-predict(v.lda)
> v.pred
> detach(v)
R程序结果：
Call:
lda(group ~ x1 + x2 + x3 + x4)
Prior probabilities of groups:
1 2 3
0.3333333 0.3333333 0.3333333
Group means:
x1 x2 x3 x4
1 188.6 150.4 13.8 20.0
2 157.0 115.0 7.0 13.6
3 151.0 121.
4 5.0 8.0
Coefficients of linear discriminants:
LD1 LD2
x1 0.01004794 0.00387986
x2 0.04017756 0.05461741
x3 0.17639800 -0.16001018
x4 0.03054815 -0.06205913
Proportion of trace:
LD1 LD2
0.9362 0.0638
$class
[1] 1 1 1 3 1 2 2 3 2 2 2 3 3 3 3
Levels: 1 2 3
结果说明：
1) Group means: 包含了每组的平均向量；
2) Coefficients of linear discriminants: 线性判别系数；
3)由$class可以看出15个训练样品判错的个数为3个（把原属于第1类的第4号观测判归为第3类；把原属于第2类的第8号观测判归为第3类；把原属于第3类的第11号观测判归为第2类）；
4) Proportion of trace: 表明了第i判别式对区分各组的贡献大小。

第六章聚类分析
聚类分析是一类将数据所研究的对象进行分类的统计方法。

这一类方法的共同点是：事先不知道类别的个数与结构，据以进行分析的数据是对象之间的相似性或相异性地数据。

将这些相似（或相异）性数据看作是对象之间的“距离”远近的一种度量，将距离近的对象归为一类，而不同类之间的对象距离较远。

6.2距离和相似系数
6.2.1距离
R中dist()函数给出了各种距离的计算结果，其使用格式为
dist(x,method="euclidean",diag=FALSE,upper=FALSE,p=2) 其中x是样本构成的数据矩阵或数据框，method表示计算距离的方法，缺省值为Euclide 距离，所定义的距离有
• "euclidean" — Euclide距离（欧式距离）；
• "maximum" — Chebyshev距离（切比雪夫距离）；
• "manhattan"—绝对值距离；
• "canberra" —Lance距离（兰氏距离）；
• "minkowski" —Minkowski距离（闵科夫斯基距离），其中参数p是Minkowski距离的阶数；
• "binary" —定性变量的距离。

diag是逻辑变量，当diag=TRUE时，给出对角线上的距离；upper是逻辑变量，当upper=TRUE时，给出上三角矩阵的值（缺省值仅给出下三角矩阵的值）。

6.2.2数据中心化与标准化变换
⑴中心化变换
⑵标准化变换
R中可用scale()函数作数据的中心化或标准化，其使用格式为：
scale(x,center=TRUE,scale=TRUE)
center=TRUE（缺省值）表示对数据作中心化变换，FALSE表示不作变换；scale=TRUE （缺省值）表示对数据作标准化变换。

⑶极差标准化变换
R中可用sweep()函数作极差标准化变换，其变换过程如下：
center <- sweep(x, 2, apply(x, 2, mean))
R <- apply(x, 2, max) - apply(x,2,min)
x_star <- sweep(center, 2, R, "/")
sweep()函数的运算格式为：
sweep(x, MARGIN, STATS, FUN="-", ...)
其中x是数组或矩阵，MARGIN是运算的区域，对于矩阵来讲，1表示行，2表示列。

STATS是统计量。

FUN表示函数的运算，缺省值为减法运算。

6.2.3相似系数
⑴夹角余弦
R中可用scale()函数完成两向量夹角余弦的计算，计算公式：
y <- scale(x, center = F, scale = T)/sqrt(nrow(x)-1)
C <- t(y) %*% y
C即是计算出的相似系数构成的矩阵。

⑵相关系数
R中利用cor( )函数即可求出样本相关阵。

6.3 系统聚类法
类与类之间的距离有许多定义方法, 主要有下面七种：。