年 月第 卷 第 期 数理统计与管理 文章编号 一 刁 刁 软件在多元统计分析中的应用 程毛林 苏州科技学院应用数学系 , 江苏苏州 , 摘要 许多实际问题往往需要对数据进行统计分析 , 建立合适的统计模型 过去一般采用 、 软件分析, 本文给出 软件在多元统计分析上的应用 , 主要介绍 在主成份分析 、 聚类分析 、 判别分析上的应用 文中均给以实例 , 结果令人满意关锐词 软件 主成份分析 聚类分析 判别分析中圈分类号 文献标识码 卜 , 盯, , , , , , , 即 滋 · 弓言 许多实际问题往往需要对数据进行多元统计分析 , 建立合适的模型 , 在多元统计分析方面 , 常用的软件有 、 、 等 我们在这里给出 在多元统计分析上的应用 , 在 较早的版本中, 统计功能不那么强大 , 而在 版本中 , 仅在统计工具箱 中的功能函数就达 多个 , 功能已足以赶超任何其他专用的统计软件 , 统计工具箱几乎包括了数理统计方面的所有概念 、 理论 、 方法 、 算法及其实现 在应用上 , 具有其他软件不可比拟的操作简单 , 接 口方便 , 扩充能力强等优势 , 再加上的应用范围广泛 , 因此可以预见 在统计应用上越来越占有极其重要的地位 下面用实例给出 在主成份分析 、 聚类分析 、 判别分析上的应用主成份分析 这里给出江苏省生态城市主成份分析实例 收稿 期 年 月 日 收到修改稿 期 年 月 日 数理统计与管理 第 卷 第 期 年 月 城市环境生态化是城市发展的必然趁势 , 表现为社会 、 经济、 环境与生态全方位的现代化水平 , 一个符合生态规律的生态城市应该是结构合理 、 功能高效和关系协调的城市生态系统所谓结构合理是指适度的人 口密度 , 合理的土地利用 , 良好的环境质量 , 充足的绿地系统 , 完善的基础设施 , 有效的自然保护 功能高效是指资源的优化配置 、 物力的经济投入 、 人力的充分发挥 、 物流的畅通有序、 信息流的快捷 关系协调是指人和自然协调 、 社会关系协调 、 城乡协调 、 资源利用和更新协调一个城市要实现生态城市的发展 目标 , 关键是在市场经济的体制 下逐步改善城市的生态环境质量 , 防止生态环境质量恶化 , 因此 , 对城市的生态环境水平调查评价很有必要 我们对江苏省十个城市的生态环境状况进行了调查 , 得到生态环境指标的指数值 , 见表现对生态环境水平分析和评价 衰 指标指数值 一级指标 结构 功能 协调 生态环境水平排二级 人口 基础 地理 城 物质 资源 生 城市 可指 结构二 设施二 结构二 绿化 还原￡。 配置二 效率劣 文明 。 续性无锡市常州镇江市 张家港 封连云港市扬州泰州徐州市南京苏州我们利用 中的 命令实现 具体程序如下二 别 妞 , · 时介 二 对武 求各变量标准差 , 。 葱 二 亡 · 亡击 , 。 , , 标准化变换 加, 那认 , 扭 二 尹讯 “几夕 调用主成分分析程序 程毛林 软件在多元统计分析中的应用 外 , 输出前三个主成分系数 二 犷讯 , 输出前三个主成分得分夕 输出特征根 。。 。。 。。 输出各个主成分贡献率执行后得到所要结果 , 这里是前三个主成分、 主成分得分、 特征根 即 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 一 夕 , , 这样 , 前三个主成分为 , , 名 一 劣 一 劣 一 。 一 一 一 劣 。 一 一 一 一 一 。 第一主成分贡献率为 , 第二主成分贡献率为 , 第三主成分贡献率

〇、前情提要
b站看【厦门大学MOOC】多元统计分析，因为老师很好看。

参考：

【厦门大学MOOC】多元统计分析

https://www.bilibili.com/video/BV1v7411E7PB

课程大纲：

一、多元统计分析概述

二、多元数据的描述与展示

三、多元正态分布

四、均值向量的检验

五、判别分析和分类分析

六、主成分分析

七、因子分析

八、聚类分析

一、多元统计分析概述


1.1 多元分析的定义
多元统计分析是什么？
多元统计分析定义

多元数据

例子
鸢尾花例子


行：样本     列：信息维度



->研究变量之间的相关性、做回归
购物网站例子




-> 维度之间的关系

-> 业务问题

1.2 多元分析的方法简介

数据描述
多元数据特征和可视化、多元正态分布

怎么从特征上面去描述：多元数据波动性、平均情况、变量与变量之间的相关性
怎么用图形表示
多元正态分布情况

第二章、第三章
统计推断
多元数据的统计检验

假设检验（数理统计 一个变量时均值的检验、两样本t检验）

第四章
经典降维
简化数据结构

具体怎么做

第六章、第七章
目标归类
根据数据特征构造归类模式

特征 聚类问题
标签 分类 判别分析问题

第八章、第五章


数据描述
将从四个部分来讲（第二章）



数据都是有分布的，多元正态（第三章）

二元正态


顾客满意度评分
平均、波动性、相关性


微博活跃程度

统计推断
μ=μ0的推广（第四章）



经典降维



不是所有的信息都有用

用少数代替多数


主成分分析
样本/个体之间差异 最大化方差
因子分析
综合指标/公共因子 变量与变量之间的相关性 有公共因素
数据减肥

目标归类
对新的样本分类



多种信息维度分类


监督学习-分类问题
分类问题



第五章


无监督学习-聚类问题
聚类问题



第八章



1.3 多元分析的应用领域
用统计学原理，研究各种感兴趣领域的知识。


不同行业的应用
聚类 分类 判别问题
市场营销

银行业

金融行业

医疗行业

分子生物学

天文学

法务会计

如何使数据驱动价值
有原始数据

->直观有效信息（二三章 可视化 波动性 平均情况 分布性）

->提取有用的知识（统计推断 显著）

->统计建模（回归 分类）

       多元统计分析是从经典统计学中发展起来的一个分支，是一种综合分析方法，它能够在多个对象和多个指标互相关联的情况下分析它们的统计规律，很适合农业科学研究的特点。主要内容包括多元正态分布及其抽样分布、多元正态总体的均值向量和协方差阵的假设检验、多元方差分析、直线回归与相关、多元线性回归与相关(Ⅰ)和(Ⅱ)、主成分分析与因子分析、判别分析与聚类分析、Shannon信息量及其应用。简称多元分析。当总体的分布是多维（多元）概率分布时，处理该总体的数理统计理论和方法。数理统计学中的一个重要的分支学科。       

       多元统计分析有狭义与广义之分，当假定总体分布是多元正态分布时，称为狭义的，否则称为广义的。近年来，狭义多元分析的许多内容已被推广到更广的分布之中，特别是推广到一种称为椭球等高分布族之中。

       按多元分析所处理的实际问题的性质分类，重要的有如下几种。

多重回归分析

简称回归分析。其特点是同时处理多个因变量。回归系数和常数的计算公式与通常的情况相仿，只是由于因变量不止一个，原来的每个回归系数在此都成为一个向量。因此，关于回归系数的检验要用T2统计量；对回归方程的显著性检验要用Λ统计量。

回归分析在地质勘探的应用中发展了一种特殊的形式，称为趋势面分析，它以各种元素的含量作为因变量，把它们对地理坐标进行回归（选用一次、二次或高次的多项式）,回归方程称为趋势面,反映了含量的趋势。残差分析是趋势面分析的重点，找出正的残差异常大的点，在这些点附近，元素的含量特别高，这就有可能形成可采的矿位。这一方法在其他领域也有应用。

判别分析

由 k个不同总体的样本来构造判别函数，利用它来决定新的未知类别的样品属于哪一类，这是判别分析所处理的问题。它在医疗诊断、天气预报、图像识别等方面有广泛的应用。例如，为了判断某人是否有心脏病，从健康的人和有心脏病的人这两个总体中分别抽取样本，对每人各测两个指标X1和X2，点绘如图 。

判别分析图

可用直线A将平面分成g1和g2两部分，落在g1的绝大部分为健康者，落在g2的绝大部分为心脏病人,利用A的垂线方向l=(l1,l2)来建立判别函数

y=l1X1+l2X2,可以求得一常数с，使 y<с 等价于（X1，X2）落在g1，y>с等价于（X1，X2）落在g2。由此得判别规则：若，l1X1+l2X2<c

判,

.

即此人为健康者;若，l1X1+l2X2>C

判，

.

即此人为心脏病人;若，l1X1+l2X2=c则为待判。此例的判别函数是线性函数,它简单方便,在实际问题中经常使用。但有时也用非线性判别函数，特别是二次判别函数。建立判别函数和判别规则有不少准则和方法，常用的有贝叶斯准则、费希尔准则、距离判别、回归方法和非参数方法等。

无论用哪一种准则或方法所建立的判别函数和判别规则，都可能产生错判，错判所占的比率用错判概率来度量。当总体间区别明显时，错判概率较小；否则错判概率较大。判别函数的选择直接影响到错判概率，故错判概率可用来比较不同方法的优劣。

变量（如上例中的X1和X2）选择的好坏是使用判别分析的最重要的问题，常用逐步判别的方法来筛选出一些确有判别作用的变量。利用序贯分析的思想又产生了序贯判别分析。例如医生在诊断时,先确定是否有病,然后确定是哪个系统有病，再确定是什么性质的病等等。

聚类分析

又称数值分类。聚类分析和判别分析的区别在于，判别分析是已知有多少类和样本来自哪一类，需要判别新抽取的样本是来自哪一类；而聚类分析则既不知有几类,也不知样本中每一个来自哪一类。例如,为了制定服装标准，对 N个成年人，测量每人的身高(x1)、胸围(x2)、肩宽(x3)、上体长(x4)、手臂长(x5)、前胸(x6)、后背(x7)、腰围(x8)、臀围(x9)、下体长(x10)等部位，要将这N个人进行分类,每一类代表一个号型；为了使用和裁剪的方便，还要对这些变量(x1,x2,…,x10)进行分类。聚类分析就是解决上述两种分类问题。

设已知N个观测值X1,X2,…,Xn,每个观测值是一个p维向量（如上例中人的身高、胸围等）。聚类分析的思想是将每个观测值Xi看成p维空间的一个点,在p维空间中引入“距离”的概念，则可按各点间距离的远近将各点（观测值）归类。若要对 p个变量（即指标）进行分类，常定义一种“相似系数”来衡量变量之间的亲密程度，按各变量之间相似系数的大小可将变量进行分类。根据实际问题的需要和变量的类型，对距离和相似系数有不同的定义方法。

按距离或相似系数分类,有下列方法。①凝聚法:它是先将每个观察值{Xi}看成一类，逐步归并，直至全部观测值并成一类为止，然后将上述并类过程画成一聚类图（或称谱系图），利用这个图可方便地得到分类。②分解法：它是先将全部观测值看成一类，然后逐步将它们分解为2类、3类、…、N类,它是凝聚法的逆过程。③动态聚类法：它是将观测值先粗糙地分类，然后按适当的目标函数和规定的程序逐步调整，直至不能再调为止。

若观察值X1,X2,…,Xn之间的次序在分类时不允许打乱，则称为有序分类。例如在地质学中将地层进行分类，只能将互相邻接的地层分成一类，不能打乱上下的次序。用于这一类问题中的重要方法是费希尔于1958年提出的最优分割法。

聚类分析也能用于预报洪水、暴雨、地震等灾害性问题，其效果比其他统计方法好。但它在理论上还很薄弱，因为它不象其他方法那样有确切的数学模型。



主成分分析

又称主分量分析，是将多个变量通过线性变换以选出较少个数重要变量的一种方法。设原来有p个变量x1,x2,…,xp,为了简化问题，选一个新变量z，

公式

,

要求z尽可能多地反映p个变量的信息，以此来选择l1,l2,…,lp,当l1,l2,…,lp选定后,称z为x1,x2,…,xp的主成分（或主分量）。有时仅一个主成分不足以代表原来的p个变量，可用q(<p)个互不相关的呈上述形式的主成分来尽可能多地反映原p个变量的信息。用来决定诸系数的原则是，在

公式

的约束下，选择l1,l2,…,lp使z的方差达到最大。

在根据样本进行主成分分析时又可分为R型分析与Q型分析。前者是用样本协差阵（或相关阵）的特征向量作为线性函数的系数来求主成分；后者是由样品之间的内积组成的内积阵来进行类似的处理，其目的是寻找出有代表性的“典型”样品，这种方法在地质结构的分析中常使用。



对应分析

这是70年代地质学家提出的方法。对非负值指标的样本资料矩阵作适当的处理后,同时进行R型与Q型的主成分分析,将结果综合在图上进行解释，可以得到指标随时间、空间位置变化的规律。它的理论正在引起多方面的重视。



因子分析

它是由样本的资料将一组变量

公式

y2，……yp）

分解为一些公共因子f与特殊因子s的线性组合,即有常数矩阵A使у=Af+s。公共因子f 的客观内容有时是明确的,如在心理研究中,根据学生的测验成绩（指标）来分析他的反应快慢、理解深浅（公共因子）；有时则是不明确的。为了寻求易于解释的公共因子，往往对因子轴进行旋转，旋转的方法有正交旋转,斜旋转,极大变差旋转等。

从样本协差阵或相关阵求公共因子的方法有广义最小二乘法、最大似然法与不加权的最小二乘法等。通常在应用中，最方便的是直接利用主成分分析所得的头几个主成分，它们往往是对各个指标影响都比较大的公共因子。



典型相关分析

它是寻求两组变量各自的线性函数中相关系数达到最大值的一对，这称为第一对典型变量，还可以求第二对，第三对，等等，这些成对的变量，彼此是不相关的。各对的相关系数称为典型相关系数。通过这些典型变量所代表的实际含意，可以找到这两组变量间的一些内在联系。典型相关分析虽然30年代已经出现，但至今未能广泛应用。

上述的各种方法可以看成广义多元分析的内容，在有些方法中，如加上正态性的假定，就可以讨论一些更深入的问题，例如线性模型中有关线性假设检验的问题，在正态的假定下，就有比较系统的结果。 　多元分析也可按指标是离散的还是连续的来区分，离散值的多元分析实质上与列联表分析有很大部分是类似的，甚至是一样的。

非数量指标数量化的理论和方法也是广义多元分析的一个重要的研究课题。

 

引言

本学期也开了一门多元统计分析课程，也趁机想把课后上机题实现一遍，以增强理解。

教材使用的是约翰逊的《多元统计分析》第六版，中英文版教材、数据集、讲义见 

还参考了王斌会老师的《多元统计分析及R语言建模》

本文内容主要为第4章多元正态分布的上机题，图略。 

[rmd文档见](http://pan.baidu.com/s/1ntkuXQT） 

可以直接用Rstudio打开（之前先安装knitr包）



4.28



data_4.28<-read.table("E:\\研究生\\应用多元统计\\JohnsonWichern Data sets\\T1-5.DAT")
#正态Q-Q图
qqnorm(data_4.28$V2)
#正态性检验
#原始数据排序
new_data<-sort(data_4.28$V2)
length(new_data)
#对应概率值

prob<-function(i,n=42){#构建一个概率值的函数
  return((i-0.5)/n)
}
all_pro<-sapply(1:42,prob)#所有概率值
#对应的标准正态分位数
all_q<-qnorm(all_pro)
#Q-Q图的相关系数
rq<-cor(new_data,all_q)
#由于Q-Q图的相关系数rq为0.9693258，小于表4-2中n=40对应的临界点，所以拒绝正态性假设。



4.29



#(a)
#计算样本协方差矩阵
s<-cov(data_4.28[,5:6])
#s的逆
s_solve<-solve(s)
x_bar<-apply(data_4.28[,5:6],MARGIN=2,mean)#两列平均数
x_bar<-matrix(as.vector(x_bar),42,2,by=2)
two_col<-t(data_4.28[,5:6]-x_bar)#两列x-x_bar
#计算所用统计距离dis
dis<-c()
for(i in 1:length(two_col[1,])){
  dis[i]<-t(two_col[,i])%*%s_solve%*%two_col[,i]
}
####################################
#(b)
#自由度为2概率密度为0.5的卡方分布临界值
chisq_num<-qchisq(0.5,2)
#所占比例
pro<-length(which(dis<chisq_num))/length(dis)
####################################
#(c)
#对广义平方距离dis进行排序
sort_data<-sort(dis)
#概率密度为4.28中的all_pro
#对应的自由度为2的卡方分位数

all_chiisq<-sapply(all_pro,qchisq,df=2)#所有概率值
#画出卡方图 也就是(all_chiisq,sort_data)对应的散点图
library(ggplot2)
qplot(all_chiisq, sort_data, geom='point')



4.30



#读入数据
data_4.30_x1<-c(1:9,11)
data_4.30_x2<-c(18.95,19.00,17.95,15.54,14.00,12.95,8.94,7.49,6.00,3.99)

#构建幂变化函数
##幂类变化函数(Box-Cox)
box_cox<-function (x,λ){
  if (λ==0) { 
    return(log(x))
  }else{
    return((x^λ-1)/λ)
  }
}
l_value<-function(X,lamda){
  x_new<-sapply(X,box_cox,λ=lamda)
  x_bar<-mean(x_new)
  l_val<-log(mean((x_new-x_bar)^2))*(-length(x_new)/2)+(lamda-1)*sum(log(X))
  return(l_val)
}

#生成多个λ，求使l_value最大的λ_hat值
λ<-seq(-1,2,0.1)
all_l<-c()
for(n in 1:length(λ)){
  all_l[n]<-l_value(data_4.30_x1,lamda=λ[n])
}
#取使变化后的l_value最大的λ值
max_λ<-λ[which(all_l==max(all_l))]
#进行数据幂变化
new_data<-sapply(data_4.30_x1,box_cox,λ=max_λ)
#变化后的Q-Q图
qqnorm(new_data)
###################################
#(b)
#基本同(a)题
λ<-seq(-1,2,0.1)
all_l<-c()
for(n in 1:length(λ)){
  all_l[n]<-l_value(data_4.30_x2,lamda=λ[n])
}
#取使变化后的l_value最大的λ值
max_λ<-λ[which(all_l==max(all_l))]
#进行数据幂变化
new_data<-sapply(data_4.30_x2,box_cox,λ=max_λ)
#变化后的Q-Q图
qqnorm(new_data)
#################################
#(c)略
#题4.31-4.38均按照4.28-4.30的解题思路进行即
#考虑边缘正态性：先做Q-Q图做个粗略的了解 然后计算Q-Q图的相关系数 并与书中表4.2进行比较 得出是否拒绝正态性的假设
#考虑二维正态性 采用4.29的方法 做卡方图
#变换可以采用平方根变换 对数变换 z变换 ，见书本p147页，还可以使用4.30中的幂变换，然后将变换后的数据画Q-Q图进行判断。



4.39



data_4.39<-read.table("E:\\研究生\\应用多元统计\\JohnsonWichern Data sets\\T4-6.DAT")[,1:5]
#(a)
#正态性检验
#计算Q-Q图的相关系数 将题4.28的代码进行封装
norm_test<-function(data){
  #原始数据排序
  new_data<-sort(data)
  len_data<-length(new_data)
  prob<-function(i,n){#构建一个概率值的函数
    return((i-0.5)/n)
  }
  #对应概率值
  all_pro<-sapply(1:len_data,prob,n=len_data)#所有概率值
  #对应的标准正态分位数
  all_q<-qnorm(all_pro)
  #Q-Q图的相关系数
  return(cor(new_data,all_q))
}
##对于独立性
#Q-Q图
qqnorm(data_4.39$V1)#大部分在一条直线上
norm_test(data_4.39$V1)
#在显著性水平为0.05的情况下，当n=150时，0.988小于于表4.2中的0.9913拒绝正态性假定。
#也可以采用shapiro-wilk检验
#使用在mvnormtest包里mshapiro.test，具体可以使用?mshapiro.test查看使用方法
##对于支撑力
qqnorm(data_4.39$V2)#大部分在一条直线上
norm_test(data_4.39$V2)
#在显著性水平为0.05的情况下，当n=150时，0.989小于表4.2中的0.9913拒绝正态性假定
##对于仁爱心
qqnorm(data_4.39$V3)#大部分在一条直线上
norm_test(data_4.39$V3)
#在显著性水平为0.05的情况下，当n=150时，0.993大于表4.2中的0.9913不拒绝正态性假定
#对于顺从性
qqnorm(data_4.39$V4)#大部分在一条直线上
norm_test(data_4.39$V4)
#在显著性水平为0.05的情况下，当n=150时，0.993大于表4.2中的0.9913 不拒绝正态性假定
#对于领导能力
qqnorm(data_4.39$V5)#大部分在一条直线上
norm_test(data_4.39$V5)
#在显著性水平为0.05的情况下，当n=150时，0.981小于表4.2中的0.9913 拒绝正态性假定
###################################
#(b)
##使用卡方图进行判定
#构造画卡方图的函数 方法同题4.29
chis_chart<-function(x){
  #计算样本协方差矩阵
  s<-cov(x)
  #s的逆
  s_solve<-solve(s)
  x_bar<-apply(x,MARGIN=2,mean)#两列平均数
  two_col<-t(x-x_bar)#两列x-x_bar
#计算所用统计距离dis
  dis<-c()
  for(i in 1:length(two_col[1,])){
    dis[i]<-t(two_col[,i])%*%s_solve%*%two_col[,i]
  }
  #对广义平方距离dis进行排序
  sort_data<-sort(dis)
  #prob在题4.28中构造
  all_pro<-sapply(1:length(x[,1]),prob,n=130)#所有概率值
  #对应的自由度为5的卡方分位数
  all_chiisq<-sapply(all_pro,qchisq,df=5)#所有概率值
  #画出卡方图 也就是(all_chiisq,sort_data)对应的散点图
  library(ggplot2)
  qplot(all_chiisq, sort_data, geom='point')
}
  chis_chart(data_4.39)
  #很明显，卡方图上点不是接近于一条直线，偏一条曲线，所以多元正态性不满足，可知，边缘正态性不满足的情况下，多元正态性也很少满足
#
###################################
#(c)
#在(a)中，独立性、支撑力、领导力的分布不符合正态性
##幂变化函数构造见题4.30
##对于独立性
#生成多个λ，求使l_value最大的λ_hat值
λ<-seq(-1,2,0.1)
all_l<-c()
for(n in 1:length(λ)){
  all_l[n]<-l_value(data_4.39$V1,lamda=λ[n])
}
#取使变化后的l_value最大的λ值
max_λ<-λ[which(all_l==max(all_l))]
#进行数据幂变化
new_data<-sapply(data_4.39$V1,box_cox,λ=max_λ)
#变化后的Q-Q图
qqnorm(new_data)
##对于支撑力
all_l<-c()
for(n in 1:length(λ)){
  all_l[n]<-l_value(data_4.39$V2,lamda=λ[n])
}
#取使变化后的l_value最大的λ值
max_λ<-λ[which(all_l==max(all_l))]
#进行数据幂变化
new_data<-sapply(data_4.39$V2,box_cox,λ=max_λ)
#变化后的Q-Q图
qqnorm(new_data)
##对于领导力
all_l<-c()
for(n in 1:length(λ)){
  all_l[n]<-l_value(data_4.39$V5,lamda=λ[n])
}
#取使变化后的l_value最大的λ值
max_λ<-λ[which(all_l==max(all_l))]
#进行数据幂变化
new_data<-sapply(data_4.39$V5,box_cox,λ=max_λ)
#变化后的Q-Q图
qqnorm(new_data)


4.40



data_4.40<-read.table("E:\\研究生\\应用多元统计\\JohnsonWichern Data sets\\T1-11.DAT")
library(ggplot2)
#散点图检查
qplot(data_4.40$V1, data_4.40$V2, geom='point')
#从散点图可以看出在x轴和y轴分别有一个离群值
#标准化值来检查
cen_data<-scale(data_4.40)
#每一列的最大离群值为
apply(abs(cen_data),2,max)
#与取标准化数据比较，第一列第13行，第二列第7行与其他数据存在较大偏离
#(b)(c)略4.40略