多元统计分析

Transcription

1 多元统计分析 ( 统计专业 ) 主讲张小霞

2 第一章概率总论 多元概率

3 0. 分布函数 离散分布 连续分布 正态分布 norm 二项分布 binom(n,p) 几何分布 geom(p) 超几何分布 hyper(m,n,k) 负二项分布 binom(x,size,prob,mu) 泊松分布 pois(lambda) 函数对应意义 d 对应概率密度函数或概率质量函数 P(X=x) p 对应累计分布函数 P(X x) q 对应分布的分位数 r 对应随机数生成函数 贝塔分布 beta(shape1,shape2, ) 柯西分布 cauchy(location,scale) 卡方分布 chisq(df) 指数分布 exp(rate) F 分布 f(df1,df2, ) 伽马分布 gamma(rate,scale) 对数正态分布 lnorm(meanlog,sdlong) 逻辑分布 logis(location,scale) 正态分布 norm(mean,sd) 学生 t 分布 t(df) 均匀分布 nuif(min,max) 威布尔分布 weibull(shap,scale) Wilcxon 分布 wilcox(m,n)

4 1.1 组合数 计算组合数 choose(n,k) 从 n 个中取 k 个的组合数 choose(60,5) # How many ways can we select 5 items from 60 items? [1] 生成组合 combn(items,k); combn(1:4,2) [,1] [,2] [,3] [,4] [,5] [,6] [1,] [2,]

5 Random Walk 随机 Time 生成随机数 runif(min=0,max=1,3) # 生成 0 至 1 之间的 3 个随机数 rnorm(mean=0,sd=1,100) # 生成标准正态分布 100 个随机数 rbinom(2,size=10,prob=0.3) # two binomial variates rpois(3,lambda=10) # 3 Poisson variates 生成可再生随机数 需要生成一个随机序列, 在程序每次运行时都复制此序列 set.seed(n) # n is a positive integer <-set.seed(15)0) <-runif(10) # 上述两个命令同时使用, 会得到相同的随机数

6 Random Walk 随机 Time 生成随机样本 sample(vec,n,replace=f) # F 为不放回抽样,T 为放回抽样 生成随机序列 sample(set,n,replace=t) sample(c(false,true),20,replace=t,prob=c(0.2,0.8)) # 生成 20 个伯努利实验, 成功概率 p= 随机排列向量 sample(v,size=length(v),replace=false) x<-1:10; sample(x,size=length(x),replace=false)

7 1.3 计算分布函数概率 对于一个累计概率,P(X x), 使用分布函数计算. 二项式 :pbinom(x,size,prob), 几 何 :pgeom(x,prob) 泊 松 :ppois(x,lambda) 正 态 :pnorm(x,mean,sd) 学生 t: pt(x,df) 指 数 :pexp(x,rate) 伽 马 :pgamma(x,shape,rate) 卡 方 :pchisq(x,df)

8 1.4 计算分布函数概率例子 二项分布 pbinom(7,size=10,prob=0.4); 7 10 数学公式为 k=0 k R 程序为 (pbinom 函数为 ) pbinom<-function(n,s,p){ x<-0 for(i in 0:n){ 0.4k k j<-i+1; x[j]<-choose(m,i)*p^i*(1-p)^(m-i)} sum(x) } pbinom(c(3,7),size=10,prob=0.4); diff(pbinom(c(3,7),size=10,prob=0.4))# 计算累计概率的差

9 1.5 绘制密度函数图形 dnorm(x) plot x<-seq(from=-3,to=3,length.out=100) plot(x,dnorm(x)) 标准正态密度函数 x x<-seq(from=0,to=6,length.out=100) ylim<-c(0,0.6) par(mfrow=c(2,2)) Plot(x,dunif(x,min=2,max=4),main= Uniform,type= 1,ylim=ylim)

10 1.5 绘制密度函数图形 1.5.1plot plot(x,dnorm(x,mean=3,sd=1),main= Normal,type= 1,ylim=ylim) plot(x,dchisq(df=5),main= Chisq,type= 1,ylim=ylim)

11 1.5 绘制密度函数图形 polygon: region 阴影区域画法 调用函数 polygon 为一个区域打上阴影 例 : x<-seq(from=-3,to=3,length=100) y<-dnorm(x) plot(x,y,main= Standard Normal Distribution,type= l,ylab= Density,xlab= Quantile ) abline(h=0) region.x<-x[1<=x&x<=2]; region.y<-y[1<=x&x<=2]

12 1.5 绘制密度函数图形 region.x<-c(region.x[1],region.x,tail(region.x,1)); region.y<-c(0,region.y,x) polygon(region.x,region.y,density=10) Standard Normal Distribution Density Quantile

13 第二章统计概论 相关 R 函数及图形

14 第三章多元正态方法 相关 R 函数及图形

15 3.1 多元分析方法概要 多元分析方法目的 多元分析方法关心的是在相互之间关系 (1) 响应变量之间 (2) 实验单元之间 (3) 响应变量和实验单元之间 Variable- and Individual- Directed Techniques 变量技术是指 : 可测量的响应变量之间主要的关系 主要集中在相关矩阵 主成分分析 因子分析 回归分析和典型相关分析 个体直接技术 : 是指实验单元之间的关系 主要集中在判别分析和典型 聚类分析和多变量分析 (MANOVA)

16 3.1 多元分析方法概要 创建新变量 我们常常发现创建新变量是非常有用的, 许多多元方法帮助研究者创建具有希望性质的新变量 例如 : 主成分分析 因子分析 典型相关分析 典型判别分析 典型变量分析 主成分分析 (PCA) 分析一个新的数据集, 以下几个问题需要考虑 : (1) 数据集是否有特殊或非一般的现象? (2) 数据是否假定为正态分布? (3) 是否有其他的非正态分布的数据? (4) 数据中是否有 outliers? 主成分分析 : 主要是创建一个不相关的数据集或随机变量称为主成分 这样的主成分是通过对向量的正交变换得来的 通过主成分得分 (principal components scores) 来判别回答上述四个问题

17 3.1 多元分析方法概要 因子分析 (FA) 因子分析技术主要是创建新的随机变量来描述原随机变量中的信息 分为公共因子和独特因子 主要研究数据集中变量之间的关系, 描绘向量之间的高相关性和低相关性 找出影响数据表象变化的内在因素, 即为公共因子 注 : 由 FA 创建的新变量公共因子比由 PCA 创建的新变量主成分好解释! 判别分析 (DA) 以例子说明 : 作为银行发放信用卡, 首先要解决把人群分为两类 :(1) 很好的信用 good credit risks(2) 信用风险高 bad credit risks 为了区分人群, 银行可能把教育水平 工资水平 债务及其信用历史作为可能的将来的信誉 (creditworthiness), 依据这些公司才能决定申请者的信用为多少, 多元统计方法能帮助公司把申请人分类的方法就叫判别分析 判别分析 : 主要是把个体或实验单元分为两类或更多类 前提为必须有随机样本建立规则 银行主要依赖于人口统计中以往记录

18 3.1 多元分析方法概要 典型判别分析 (DFA) 是判别分析创建有用信息的新变量预处理 新变量对分不同类有着简单的规则 Logistic Regression 是一个概率模型, 这个模型主要应用于判断 在信用卡的例子中, 可以模拟这样的模型做出判断 聚类分析 Cluster analysis 聚类分析 CA 类似于判别分析来分类 当研究者预先通过随机样本有了一定的子类后用判别分析 ; 而聚类分析在不知道任何信息时所用的分析方法

19 3.1 多元分析方法概要 多元方差分析 (MANOVA) 推广了一元的方差分析 主要技术为在测量一个随机变量时, 比较它们的均值 典型相关分析 Canonical Correlation analysis

20 3.2 多元正态分布 The Multivariate Normal Distribution 定义 x 1 一个随机向量 x = x 2 x p 如果满足 a x = a 1 a 2 a p x 1 p x 2 = i=1 a i x i x p 对任何不全为 0 的 a i 都是单元正态分布, 我们称向量 x 为多元正态分布 注 : 如果 x i 都是正态分布, x 不一定就 Mean vectors 和 Variance- Covariance Matrices Mean E(x 1 ) μ = E(x 2 ) E(x p ) μ 1 = μ 2 μ p 协变矩阵 =Cov(x)=E[(x-μ)(x- μ) ] = σ 11 σ 12 σ 1p σ 21 σ 22 σ 2p σ p1 σ p2 σ pp σ ii = Var x i = E[(x i μ i ) 2 ] σ ij = Cov x i, x j = E[ x i μ i x j μ j ] 是多元正态分布 反例见 Hogg and Craig(1978,p.114, Ex. 3.55)Intruoduction to Mathematical Statistics(4 th ed.) New York: Macmillan,1978

21 3.2 多元正态分布 The Multivariate Normal Distribution x 1 与 x 2 相关系数定义 ρ ij = σ ij σ ii σ jj 相关系数阵 P 1 ρ 12 ρ 21 1 ρ 1p ρ 2p ρ p1 ρ p2 1 如果 ρ ij 接近于 1, 我们称它们线性相关或简称相关 如果 ρ ij 接近于 0, 我们称它们线性无关 uncorrelated( 等价于它们组成的联合正态分布是独立的 ) 或简称无关 注 : 指的无关 uncorrelated 并不是我们所说的一点关系没有 unrelated. x_1 为均值为 0, 方差为 8 的正态分布 x_2=x_1^2-2*x_1+3 σ 12 =Cov(x_1,x_2)=0

22 3.2 多元正态分布 The Multivariate Normal Distribution 多维正态分布密度函数 x~n p (μ, ) 如果 的秩为 p, 则 x 的密度函数存在, 并密度函数为 f x x; μ, = 1 (2π) p 2 对任意的 x E p 1 2 exp[ 1 2 x μ 1 (x μ)] 其中 E p = x: < x i <, for i = 1,, p 当 p=2 时 ; 右边为其图像 画三维图像有 3 种方法 scatterplot3d(x,y,z) 点 线 persp(x,y,z) 曲面 rgl packages 动态

23 3.3 多元统计量 --- 估计量 数据统计中的一些估计量 设 x 1, x 2,, x N 是从有均值向量 μ 协方差阵 和相关矩阵 P 的多元分布中的 N 个样本 μ 和 的无偏估计为 (1) μ = 1 N (2) = 其中 μ = 1 N 1 x 1 x p = i=1 N x i = x 1+x 2 + +x N N i=1 N (x i μ)(x i μ) μ 1 μ p 相关系数 ρ ij 的估计用 r ij = σ ij σ ii σ jj 样本相关阵为 1 r 12 r 13 r 1p r 21 1 r 23 r 2p R = r 31 r 32 1 r 3p r p1 r p2 r p3 1 = σ 11 σ 1p σ p1 σ pp σ ii = σ ij = 1 N 1 1 N 1 r=1 N (x ri x i ) 2 r=1 N (x ri x i )(x rj x j ), i j

24 3.4 标准数据 and/or scores 为了比较数据和理解数据, 常常把数据标准化 最常用的就是标准化变换 z rj = x rj μ j σ jj The variable z rj is called the Z score for the jth response variable on the rth experimental unit, and Z = z z 11 z 12 z 13 1p z 21 z 22 z 23 z 2p z 31 z 32 z 33 z 3p z p1 z p2 z p3 称之为 Z 的矩阵 z pp 注 : 当数据单位不一致时, 建议标准化 R 程序 ( 从数据 M 开始 )##R 函数 core score<-function(m){ p<-ncol(m) n<-nrow(m) B<-0 mu<-0 for(i in 1:p){ mu[i]<-mean(m[,i]) } for(i in 1:n){ B<-B+1/(n-1)*(M[i,]- mu)%*%t(m[i,]-mu) } Z<-matrix(0,ncol(B),ncol(B)) for(r in 1:ncol(B)){ for(i in 1:ncol(B)){ Z[r,i]<-(M[r,i]-mu[i])/sqrt(B[i,i]) } } Z }

25 第三章例题 1 考察数据矩阵 x = 回答问题 a. p N x 32 和 x 3 的值? x 3? b. 计算均值 μ, 样本协方差阵和样本相关系数阵 R c. 计算 score 矩阵 Z d. z 32 z 3 的值? e. 画出变量 x 3 和 x 2 的散点图 f. 画出变量 z 3 和 z 2 的散点图 要求 : 以上的 b c e f 用 R 语言, 答案提示 : μ 相应 R 程序如下 mu<-function(m){ n<-ncol(m) mu<-0 for(i in 1:n){ mu[i]<-mean(m[,i]) } mu }

26 第三章例题提示 计算协方差阵 程序 sigma<-function(m){ n<-nrow(m) m<-ncol(m) B<-matrix(0,m,m) mu<-0 for(j in 1:m){ mu[j]<-mean(m[,j])# 求均值向量 } for(i in 1:n){ B<-B+1/(n-1)*(M[i,]-mu)%*%t(M[i,]-mu) } B } 计算相关系数阵 R 程序 #B 为协方差阵 R<-function(B){ R<-matrix(0,ncol(B),ncol(B)) for(i in 1:ncol(B)){ for(j in 1:ncol(B)){ R[i,j]<-B[i,j]/sqrt(B[i,i]*B[j,j]) } } R }

27 第四章样本相关性 Sample Correlations 建立在系数阵上的随机变量的分类

28 综 述 从大量的实验中有很多的不同类型的数据, 当然我们要问不同类型数据之间是否相关? 它们的相关性是怎样的关系? 为考察这样的相关性, 我们一般先把数据进行分类, 不同类之间的变量有比较低的相关性 这样就有可能揭示数据之间的现象 考察下面的例子 48 个人申请一个大公司的一个 offer, 它们进过面试及对设计的 15 项目进行了打分 15 项为 : 1. 申请信格式 (FL)2. 外貌 (APP), 3. 学术能力 (AA),4. 魅力 (LA), 5. 自信力 (SC), 6. 外向度 (LC), 7. 诚实度 (HON), 8. 说服力 (SMS), 9. 经历 (EXP), 10. 驾驶水平 (DRV) 11. 追求 (AMB), 12. 理解力 (GSP), 13. 潜力 (POT), 14. 团队 (KJ), 15. 适应性 (SUIT). 每一项都有一个从 0 到 10 的分值可填 0 指非常的低及非常不满意, 10 很高的评价 假设公司提供 offers 选取 6 个最好的申请者? 公司应该怎样选取这 6 个人? 数据见 MSA---offer 数据

29 offer 方法的讨论 1. 最简单的方法, 是不管上述的打分, 随机选取, 人人有意见! 2. 求平均 : 大家都认为每一项都同等重要!? 这是不可能, 只是为了简单办事 3. 加权平均, 这是大家所希望的 μ = 1 15 ω 1FL + ω 2 APP + + ω 15 SUIT, ω 1 + ω ω 15 = 1 在加权平均中系数怎样确定? 通过 R 语言中函数或我们编写的函数计算相关系数阵为如下图表示 ( 学生上机计算 ) 取两位小数

30 offer 相关系数阵 FL APP AA LA SC LC HON SMS EXP DRV AMB GSP POT KJ SUIT FL APP AA LA SC LC HON SMS EXP DRV AMB GSP POT KJ SUIT

31 data<-read.table( 路径 ) # 读取文件 M<-cor(data) # 变为矩阵 写一个判断列中协方差大于等于 0.8, 小于 1 的函数 f<-function(x,k1,k2){ Which(x<k2 & x>=k1) } apply(m,2,f,k1,k2)

32 APP AA LA HON POT KJ LC SC SMS DRV AMB GSP POT SC LC SMS DRV AMB GSP POT HON LA EXP SUIT DRV SC LC SMS AMB GSP POT KJ SUIT AMB SC LC SMS DRV GSP POT GSP SC LC SMS DRV AMB POT POT LA SC LC SMS DRV AMB GSP KJ LA DRV

33 具体分组规则 ( 讨论是否合理 ) 分组规则 : 相关性比较大的在一组 ( 是否有主观因素?) 第一组 :LC GSP POT SC SMS DRV AMB 第二组 :FL EXP SUIT 第三组 :LA HON KJ 第四组 :AA 第五组 :APP 问题 (1) 编程序自动的分组, 然后人工的调整? (2) 我们这种分组对求职的人是否有影响? 决策 : 在开始时, 设计的是 15 项特征, 当前的分组是设计 5 个方面不同的特征 求职者应在在这 5 个不同的特征上进行赋值, 来决定是否给予 offer. 所以在 5 方面的总分或者为平均分或加权平均 这样可能更合理一些 从新加权平均为

34 最后方案 F 1 = (SC+LC+SMS+DRV+AMB+GSP+POT) 7 F 2 = FL+EXP+SUIT 3 F 3 = LA+HON+KJ 3 F 4 = AA F 5 = APP scores = ω 1 F 1 + ω 2 F 2 + ω 3 F 3 + ω 4 F 4 + ω 5 F 5 注 : 上面的例子的处理方法是因子分析的一个很粗略地形式, 由以前的 15 个特征变为新的 uncorrelated 的 5 个特征, 这些新特征是求职者潜在的特征

35 第 5 章主成分分析 Principal components analysis First Step 数据分析师常用的经验法则 : 主成分和因子分析需要 5~10 倍于变量数的样本数

36 5.1 美国海军收集了一组数据, 目的想对学士学位毕业的学生估计人力资源的需求 数据间见 navy_data.txt Myers(1990) 用上述数据建立了多元回归, 说明了 MMH 与前面变量的关系 在进行多元回归之前, 要判断下面三个问题 (1) 数据集中是否有离群点? 因为一个离群点, 可能很戏剧的影响回归分析! (2) 前面的各变量中是否存在多重线性关系? 因为多重线性关系的存在对合适的回归模型的解释有很大的影响! (3) 考察的这些变量是否为多元正态分布? 如果不是, 我们通过变化, 变为多元正态分布 主成分分析可以对上述问题给出答案!

37 5.1 应用主成分分析的理由 1. 数据的筛选 : 建议 PCA 是作为数据分析的第一步, 通过 PCA 可以发现离群点 2. 聚类的帮助 : 把大的数据集分为小数据集合 (clustering method 涉及 ) 3. 判别分析 : 因为判别分析需要大量样品, 但遗憾的是, 往往有少量数据样品, 此时主成分作为新变量可以作为判别分析的数据样品, 从而增加数据样品的数量 4. 回归分析 : 我们知道, 当变量之间高度相关时, 多元的回归分析是非常 dangerous. 主要是中间存在多重线性 PCA 能够判别变量之间是否存在多重线性

38 5.3 基于方差 --- 协变矩阵的主成分分析 主成分性质要求 1. 所有主成分是 uncorrelated 符号意义 2. 第一主成分是所有主成分中方差最大的 3. 以后的次之 由于在数据分析中数据的均值及数据的协方差阵都是可以得到的, 而这两个多事随机变量的两个无偏差估计, 所以基于平均及协方差矩阵是已知的作为条件 这样做的前提是变量为正态分布 x 是 p 维随机向量, x 1 x = 其中 x i,i = 1, 2,, p 为随机变量 x p μ 为 x 的数学期望向量 主成分表达式 : Y 1 = a 11 X 1 + a 12 X a 1p X p Y 2 = a 21 X 1 + a 22 X a 2p X p Y p = a p1 X 1 + a p2 X a pp X p 其中 X i 为可观测的属性, 可理解为 p 个随机变量

39 为便于计算我们作下面的变化, 设 μ i 为随机变量 X i 的期 望,i = 1,, p Y 1 = a 11 (X 1 μ 1 ) + a 12 (X 2 μ 2 ) + + a 1p (X p μ p ) Y 2 = a 21 (X 1 μ 1 ) + a 22 (X 2 μ 2 ) + + a 2p (X p μ p ) Y p = a 11 (X 1 μ 1 ) + a 12 (X 2 μ 2 ) + + a 1p (X p μ p ) 简写为 Y = A X μ, 即 Y i = a i X μ, i = 1,, p

40 5.3 基于方差 --- 协变矩阵的主成分分析 (1) 理论上的简单推导 : 第一主成分的获得 : 作变换 y 1 = a 1 (x μ), 使得方差最大, cov y 1, y 1 = cov a 1 x μ, a 1 x μ = a 1 cov x μ, x μ a 1 = a 1 a 1 最大 λ 由二次型理论可知存在正交矩阵 P 使得 a 0 λ 1 P Pa 1 = b b λ p 在条件 = b b b 2 p = 1 下, 最大值为特征值最大者 λ 1 a 1 可以取协方差阵最大特征根对应的特征向量

41 第二主成分的获得 : (1) 由上述证明过程可知存在正交矩阵 P 使得 P P = λ λ λ p 若 P 写成列向量表达形式 P = a 1 a 2 a p

42 对上式变形为 λ a 1 a 2 a p = a 1 a 2 a 0 λ p λ p a 1 a 2 a p = λ 1 a 1 λ 2 a 2 λ p a p 即 a i = λ i a i, i = 1,2,, p. 结论 : 正交矩阵的列为对应的特征向量

43 5.3 基于方差 --- 协变矩阵的主成分分析 (2) (2) 条件分析 : 条件 y 1 与要求的 y 2 是 uncorrelated, 即 cov y 1, y 2 = 0 如果记 y 2 = b (x μ), 则 λ cov y 1, y 2 = a 1 b = a 1 a 1 a 2 a 0 λ p λ p a 1 a 2 a p b λ = a 1 a 1 a 1 a 2 a 0 λ 1 a p λ p a 1 b a 2 b a p b 由于 a 1 与 a i (i 1) 正交, 所以 a 1 a i = 0 (i 1), 上式变为 a 1 a λ λ λ p a 1 b a 2 b a p b = a 1 a 1 λ 1 a 1 b = 0, 从而 a 1 b = 0.

44 5.3 基于方差 --- 协变矩阵的主成分分析 (3) (3) 计算第二主成分, 假设条件 y 1 与要求的 y 2 是 uncorrelated, λ cov y 2, y 2 = b b = b a 1 a 2 a 0 λ p λ p λ λ = b a 1 b a 2 b a p λ p a 1 b a 2 b a p b = 0 b a 2 b a p 0 λ λ p λ = 0 b a 2 b a p 0 λ = b a 2 b a p 0 0 λ λ p 求上式的最大值, 就化为了类似求第一主成分的方法 0 0 λ p a 2 b a p b a 1 a 2 a p 利用 b a 1 = 0 0 a 2 b a p b 0 a 2 b a p b b

45 5.3 基于方差 --- 协变矩阵的主成分分析 (3) 主成分得分 为统计分析, 应用主成分变量 ( 把主成分作为坐标轴 ), 来计算每一个实验单元的主成分得分, 这些数据提供了对数据集的局部观察 设 x r 为第 r 个观测向量, 也可以说第 r 个实验数据项, 那么第 j 个主成分对其的得分为 y j r = a j (x r μ), for j = 1,2,, p and r = 1,2,, N. 用表格表示出来的程序如下

46 5.4 主成分估计 ---R 命令 (1) 分析流程 加载 R 包 :psych 常见步骤如下 (1) 数据处理 :PCA 和 EFA( 探索性因子分析 ) 是根据观测变量间的相关性来推导结果 用户可以输入原始数据或相关系数阵到函数 principal() 如输入初始数据, 相关系数矩阵将会自动计算, 在计算前确保数据中没有缺失值 (2) 选择模型 判断 PCA( 数据降维 ) 还是 EFA( 发现潜在结构 ) 更符合你的研究目的 (3) 判断要选择的主成分 / 因子数目 (4) 选择主成分 / 因子 (5) 旋转主成分 / 因子 (6) 解释结果 (7) 计算主成分或因子得分 我们估计中常用的数据为求的特征值和特征向量代替的特征值和特征向量

47 主成分分析步骤及相关 R 函数 1 数据分标准化 2 求相关系数矩阵 R 3 求相关系数矩阵 R 的特征根及相应的特征向量 4 按主成分贡献率确定主成分的个数, 并写出主成分表达式 5 对分词的结果做统计意义和实际意义解释

48 R 函数 1 princomp: princomp(x,cor=false,scores=true) x 为数据框或数据矩阵 ;cor=t 表示用相关系数阵, 否则用协方差阵 ;scores 是否输出主成分得分 2 summary:summary(object,loadings=f) 用于提取主成分信息, 其中 object 是指上述对象 ; loadings=t 表示显示载荷内容 3 loadings(object): 显示载荷内容, 也就是正交矩阵内容 4 predict(object): 预测主成分分析的值 5 screeplot(object,type= line,barplot )

49 综合案例 某市 13 个行业 8 项重要经济指标数据, 其中 X1 为年末固定置产净值 ( 单位 : 万元 );X2 为职工 人数数据 ;X3 为工业总产值 ( 单位 : 万元 ) #X4 为全员劳动生产率 ( 单位 : 元 / 人年 );X5 为百元固定资产原值实现产值 ( 单位 : 元 );X6 为资金利税率 (%); #X7 为标准燃料消费量 ( 单位 : 吨 );X8 为能源利用效率 ( 单位 : 万元 / 吨 ); 进行主成分分 析 name X1 X2 X3 X4 X5 X6 X7 X8 冶金 电力 煤炭 化学 机器 建材 森工 食品 纺织 缝纫 皮革 造纸 文教

50 name<-industry[1]# 查看第一列 M<-industry[-1] # 去掉第一列 M<-as.data.frame(M)# 作为数据框 industry_object<-princomp(m,cor=t,scores=t)# 作主成分分析 summary(industry_object,loadings=t) loadings(industry_object) scores<-industry_object$scores mode(scores) scores<-as.matrix(scores) rownames(scores)<-name round(scores,digits=3)

51 主成分得分 Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 冶金 电力 煤炭 化学 机器 建材 森工 食品 纺织 缝纫 皮革 造纸 文教

52 因子分析 1904 年 Charles Spearman 发表著名论文 对智力测验得分进行统计分析 标志着因子分析的起点 因子分析由于计算机的高速发展在经济学 社会学 生物学 医学 心理学及教育学等方面取得很大的成功

53 问题提出 1 为了解学生的学习能力, 观测了 n 个学生的 p 个科目的考试成绩 用 X 1,, X p 表示 p 个科目 ( 例如 : 代数 分析 几何 语文..), X (i) = x i1,, x ip, i = 1,2,, n, 表示第 i 个学生的 p 科的成绩, 根据这些数据, 哪些因素决定学生的学习能力 简单的假设, 我们假设学生学习的能力是一个不可测的随机变量, 它对成绩起着一定的作用, 可以用下面的式子, 各科目观察学生的成绩来说是一个随机变量, 从而各科目有两部分组成, 一部分为学习能力 f, 另一部分为特殊影响 ( 特殊因子 )ε i, 公式如下 X i = a i f + ε i, i = 1,, p 我们称 f 为公共因子,ε i 为特殊因子 这是最简单的因子模型

54 问题提出 2 某公司对 100 名招聘人员的知识和能力进行测试, 出了 50 道题的试卷, 其内容包括的面较广, 但总的来说归纳为 6 个方面的能力, 语言表达能力 逻辑思维能力 判断事物的敏捷和果断程度能力 思想修养 兴趣爱好 生活常识等 我们将每一个方面称为因子 显然这里所说的因子不同于回归分析中的因素, 回归分析中有比较明显的实际意义, 如人口密度 工业总产值 产量等 假设 100 人测试的分数用 X i, i = 1,2,, 100, 可以用上述因子表示成线性函数, X i = a i1 f 1 + a i2 f a i6 f 6 + ε i, i = 1,2,, 100 其中 F 1, F 2,, F 6 表示 6 个上述因子, 它们对每个 X i 都有贡献, 称为公共因子 (common factor), 它们的系数 a i1,, a i6 称为因子载荷, 它们表示应试人员在六个因子方面的能力 ε i, i = 1,.., 100, 表示应试人员的能力和知识不能被前面六个因子包括部分

55 问题总结 通过以上两个模型, 我们看到, 因子模型主要应用两个方面, 一是寻求基本结构, 简化观测系统, 将错综复杂关系对象 ( 变量或样品 ) 综合为少数几个因子 ( 不可观测的随机变量 ), 以再现因子与原始变量的内在关系 ( 类似主成分分析 ); 二是用于分类, 对 p 个变量或样品进行分类 上述例子也说明了 : 因子分析分为 R 型 ( 对随机变量作因子分析 ) 和 Q 型 ( 对样品作因子分析 )

56 预备知识 : 对称矩阵的谱分解 对任何的对称矩阵存在正交矩阵 P 使得 P P = λ λ = P λ p λ λ P λ p 若 P 写成列向量表达形式 P = a 1 a 2 a p

57 = α 1 α α 2 p 0 λ λ λ p α 1 α 2 α p =λ 1 α 1 α λ p α p α p

58 6.2 R 型 --- 因子模型 假设已知条件为可观测的随机向量 X = X 1, X 2,, X p, E X = μ, D X = Cov X μ, X μ = 我们不可观测的或需要考察的公共因子向量为 F = F 1, F 2,, F m, 公共因子相互统计无关, 即 Cov F i, F j = 0, i j, E F = 0 特殊因子向量为 ε = (ε 1, ε 2,, ε p ), 假设特殊因子之间, 特殊因子与公共因子之间统计无关, 即 Cov ε i, ε j = 0, i j, Cov F i, ε j = 0 E ε = 0, 则有 D ε = diag(σ 1 σ 2 σ p )=D X i μ i = a i1 f 1 + a i2 f a im f m + ε i, i = 1,2,, p (1) 上述模型称为正交因子模型

59 正交因子模型要解决的问题 技术层面 : (1) 确定公共因子之前的加权, 即因子载荷 a ij, i = 1,2,, p, j = 1,2,, m (2) 确定特殊因子的方差 σ i, i = 1,2,, p 应用层面 : 解释公共因子的实际意义及可观测随机变量在公共因子的体现

60 正交因子技术问题 正交因子模型的矩阵形式 X = μ + AF + ε 其中 A = (a ij ) p m 称为因子载荷矩阵 由正交因子的假设条件, 我们有 D F = Cov F, F = I m

61 正交因子技术问题 演算 考察随机向量 X 的协方差阵 = Cov(X, X) = Cov(μ + AF + ε, μ + AF + ε) = Cov(AF, AF) + Cov(ε, ε)+cov(af, ε) + Cov(ε, AF) # 线性 = ACov(F, F)A + Cov(ε, ε)+acov(f, ε) + Cov ε, AF A # 矩阵提到协变符号外 = AA + D # 条件

62 正交因子技术问题 结论 = AA + D 我们要找的参数满足上述方程 其中左边是已知的, 右边未知

63 公共因子与可观测随机变量关系 Cov X, F = Cov μ + AF + ε, F = Cov AF, F + Cov ε, F = A

64 因子模型的协方差结构 Cov X, F = AA + D = A

65 6.3 因子模型的估计 要建立正交因子模型, 就要估计载 荷矩阵 A 及特殊的方差 σ i i = 1,, p, 常用的估计方法有 : 主成分法 极大似然法和主因子法

66 6.3.1 主成分法算法 根据对称矩阵的谱分解 = λ 1 α 1 α λ p α p α p 再根据贡献率的要求确定 m λ 1 α 1 λ 1 α 1, λ 2 α 2,, λ m α m λ 2 α 2 λ m α m + diag( AA )

67 主成分法模型参数估计 A = λ 1 α 1, λ 2 α 2,, λ m α m p m D = diag( -AA ) m 的确定 : m i=1 p λ i / i=1 λ i 80%

68 主成分法 R 程序实现 offer<-read.table("c:\\users\\dell\\desktop\\offer.txt") S<-cor(offer) factor.analy1<-function(s,m){ p<-nrow(s) eig<-eigen(s) A<-((eig$vectors)[,1:m])%*%diag(sqrt(eig$values)[1:m]) D<-1-diag(A%*%t(A)) list(loadings=a,specific=d,,method=c("principle Component Method")) } factor.analy1(s,3)?? 加上名字? 因子及随机变量名 (paste(factor,1:m,sep=""))

69 $loadings [,1] [,2] [,3] [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [10,] [11,] [12,] [13,] [14,] [15,] $specific [1] [8] [15] $method [1] "Principle Component Method"

70 6.3.2 主因子法算法 主因子法是对主成分法的修正 对标准数据而言 R = AA + D R R D = AA R 称为简约相关阵 (reduced correlation matrix)

71 算法如下 : 从特殊协方差阵出发 (1) 给 D 一个合适的初始值 (σ 1 1,, σ 1 p ) (2) 按主成分法对 R 进行谱分解 λ 1 α 1 R = λ 1 α 1, λ 2 α 2,, λ p α p λ 2 α 2 λ p α p (3) 则下述的 A 1 与 D 1 为一个近似解 A 1 = λ 1 α 1, λ 2 α 2,, λ m α m D 1 = diag(i AA ) 即 :σ i = 1 h i 2 = 1 1 m a ij 2. h i 2 称为随机变量 X i 的共同度 重复以上过程

72 重复第二过程 (1) 从 D 1 出发. (2) 在根据 R = R D 1 = AA, 进行谱分解 (3) 求近似解 D 2 与 A 2 一直进行下去, 满足条件 : 循环次数或 D n D n+1 < ε, 停止

73 特殊协方差阵初始值常用取法 (1) σ i 2 = 1/r ii, 其中 r ii 为 R 1 的第 i 个对角线元素 (2) 取 h i 2 = max j i r ij, σ i 2 = 1 h i 2. (3) 当 h i 2 =1, 取 σ i 2 = 0.

74 程序如下 offer<-read.table("c:\\users\\dell\\desktop\\offer.txt") R<-cor(offer) invr<-solve(r) # 求 R 的逆 D<-1/diag(invR)

75 k<-0 repeat{ R1<-R-D;D1<-D eig<-eigen(r1) val<-eig$values vec<-eig$vectors A<-(vec[,1:3])%*%(diag(sqrt(val[1:3]))) D<-1-diag(A%*%t(A)) D<-diag(D) k<-k+1 if(sqrt(sum(d-d1)^2)<1e-5 k>20) break }

76 因子载荷 A finala<-round(a,digits=3) [,1] [,2] [,3] [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [10,] [11,] [12,] [13,] [14,] [15,]

77 整个过程写成函数 factor.analy2<-function(r,m,d){ repeat{ R1<-R-D;D1<-D k<-0 eig<-eigen(r1) val<-eig$values vec<-eig$vectors A<-(vec[,1:3])%*%(diag(sqrt(val[1:3]))) D<-1-diag(A%*%t(A)) D<-diag(D) k<-k+1 if(sqrt(sum(d-d1)^2)<1e-5 k>20) break } List(loadings=A,espefic=D,Method= Principal Factor Methed ) }

78 6.3.3 因子正交旋转 满足因子方差结构方程 = AA + D 载荷矩阵不唯一 A = AO,O 为正交矩阵 方差最大的正交旋转 使用的 R 函数为 varimax(x,eps=1e-5) X 为因子载荷矩阵,eps 为迭代终止精度

79 Offer 例子说明 R<-cor(offer) invr<-solve(r) # 求 R 的逆 D<-diag(1/diag(invR)) k<-0 repeat{ R1<-R-D D1<-D eig<-eigen(r1) val<-eig$values vec<-eig$vectors A<-(vec[,1:3])%*%(diag(sqrt(val[1:3]))) D<-1-diag(A%*%t(A)) D<-diag(D) k<-k+1 if(sqrt(sum(d-d1)^2)<1e-5 k>20) break } A A.v<-varimax(A,eps=1e-5)## 旋转 A.v

80 结果比较 旋转后 [,1] [,2] [,3] [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [10,] [11,] [12,] [13,] [14,] 旋转前 [,1] [,2] [,3] [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [10,] [11,] [12,] [13,] [14,]

81 6.4 因子得分 因子得分有两种方法 : 加权最小二乘法 回归方法 到目前为止, 计算因子得分的两种估计方法到底哪个好? 还没有结论 X = μ + AF + ε; 对每个观测给出因子估计 F 公共因子的估计即为计算因子得分

82 Bartlett s Method 或称加权最小二乘法 首先建立模型后, 即求出 A 与 D 为找到 f,bartlett 建议使得下面 在 X i 均值为 0 的前提下 X i Af D 1 (X i Af) 最小, 其中 X i 为第 i 个观测向量 解为 f i = (A D 1 A) 1 A D 1 X i, i = 1,2,, n. f i = A D 1 A 1 A D 1 (X i μ)

83 Thompson 回归法回归法 X X = Af + ε 解为 f = (I + A D 1 A) 1 A D 1 X X f i = A D 1 A 1 A D 1 (X i μ)

84 R 函数 factanal(x, factors, covmat = NULL, scores = c("none", "regression", "Bartlett"), rotation = "varimax ") X: 数据矩阵或数据框 factors: 因子个数 scores : 得分采取方法 rotation: 是否旋转

85 offer<read.table("c:\\users\\dell\\desktop\\offer.txt") offer.f<factanal(offer,factors=3,rotation= varimax,scor es= regression ) offer.f$scores

86 7 判别分析 7.1 距离判别法 7.2 Fisher 判别法 7.3 Bayes 判别法

87 7.1 距离判别法 判别分析在知道样品分类的前提下, 给新的样本按某种分类准则进行归类 距离判别法 : 理解非常简单, 距离决定归属 我们常用的距离为欧氏距离和马氏距离 设 x = (x 1, x 2,, x p ) 和 y = (y 1, y 2,, y p ) 是两个随机向量, 它们有相同的协方差矩阵 d x, y = (x y) 1 (x y) 马氏距离 mahalanobis(x,center,cov) x: vector or matrix of data with, say, p columns. cov: 协方差阵

88 7.1.1 两个总体的距离判别 给定两个总体 G 1 和 G 2, 均值分别为 μ 1 μ 2, 有相同方差, 给定一个样品, 判断它属于哪一个总体? d x, μ 1 > d x, μ 2, 则 x G 2. 简单用马氏距离演算 : d x, μ 2 1 d x, μ 2 2 = 2(x μ 1+μ 2 ) 1 (μ 2 2 μ 1 ) W x = (x μ 1+μ 2 2 W(x)>0, x G 2, W(x)<0, x G 1. ) 1 (μ 2 μ 1 ) = (μ 2 μ 1 ) 1 (x μ 1+μ 2 ) 2

89 两个总体的实际操作 X (1) = 1 n 1 X i (1) ; X (2) = 1 n 2 X i (1) 协方差矩阵的联合无偏估计为 : = 1 n 1 +n 2 2 (S 1 + S 2 ) S 1 = (X i 1 X (1) ) (X i 1 X 1 ) 当两个总体的均值有显著差异时才可以

90 7.1.2 多个个总体的距离判别 给定多个总体 G i, 均值分别为 μ i, 有相同方差, 给定一个样品, 判断它属于哪一个总体? 判断哪个最小?

91 多个总体的实际操作 X (k) = 1 n 1 X i (k), k = 1,, p 协方差矩阵的联合无偏估计为 : = 1 n 1 +n 2 + +n p p (S 1 + S S p ) S k = (X i k X (k) ) (X i k X k )

92 ir<-iris[,-5] n1<-nrow(ir[1:50,]) n2<-nrow(ir[51:100,]) n3<-nrow(ir[101:150,]) ## 均值 mu1<-apply(ir[1:50,],2,mean) mu2<-apply(ir[51:100,],2,mean) mu3<-apply(ir[101:150,],2,mean) ## 协方差阵 S1<-(n1-1)*cov(ir[1:50,]) S2<-(n2-1)*cov(ir[51:100,]) S3<-(n3-1)*cov(ir[101:150,]) ## 联合协方差的无偏差估计 Sigma<-(S1+S2+S3)/(n1+n2+n3-3) M<-matrix(0,3,150) M[1,]<-A M[2,]<-B M[3,]<-C class<-apply(m,2,which.min) A<-mahalanobis(ir,mu1,Sigma) B<-mahalanobis(ir,mu2,Sigma) C<-mahalanobis(ir,mu3,Sigma)

93 当 y = 1 (μ 1 μ 2 )x 时 (μ 1y μ 2y ) 2 Var(y) = [at (μ 1y μ 2y )] 2 a T 最大 a

94 7.2 Fisher 判别法 Fisher 于 1936 年提出该方法, 这是判别分析奠基性工作 该方法的主要思想是通过将多维数据投影到一维直线上, 使得同一类别中的数据在该直线上尽量靠拢, 不同种类数据尽量分开 从方差分析角度就是 : 组内变差尽量的小, 组间变差尽量的大 然后利用前面的距离判别法来建立判别标准 在二维平面上有两类点, 小圆圈为 G, 大圆圈 F, 按原来的坐标 x1 和 x2 很难将它们区分开, 但若把它们投影到直线 y 上, 则它们 的投影明显分为两组, 同类聚在一起, 容易区分 y 为 x1 与 x2 的线性组合 y = a 1 x 1 + a 2 x 2. 作为 p 维空间 y = a T x, 其中 a 为 p 维向量 y 相对 G1 与 G2 均值为 μ 1y = E y x G 1 = a T μ 1 ; μ 2y = E y x G 2 = a T μ 2, 方差 var(y)=a T a μ 1,μ 2 为两类 G 1,G 2 总体的均值, 它们有相同的协方差阵

95 从刚才的分析, 使得相对均值 μ 1y 与 μ 2y 的距离越大线性组合越好 ( 但这里没有体现内紧 外离, 只是从均值向量来描述 ) 利用均值向量描述, 考察比值 (μ 1y μ 2y ) 2 var(y) = [at (μ 1 μ 2 )] 2, a T a 问题简化为 : 如何选取 a, 使得上式的值最大 可以证明 : 当 a = c 1 (μ 1 μ 2 ) 时, 上式最大, 当 c = 1 时, 线性函数 y = a T x = (μ 1 μ 2 ) T 1 x 为 Fisher 线性判别函数

96 相对在 y 直线上的两个均值 μ 1y, μ 2y, 它们的中间值为 μ y = 1 2 μ 1y + μ 2y = 1 2 (μ 1 μ 2 ) T 1 (μ 1 + μ 2 ) μ 2y μ y μ 1y Fisher 判别准则 : y = (μ 1 μ 2 ) T 1 x μ y, 则 x G 1 y = (μ 1 μ 2 ) T 1 x < μ y, 则 x G 2 μ 1 μ 2 T 1 (2x μ 1 μ 2 ) 0, 则 x G 1 μ 1 μ 2 T 1 (2x μ 1 μ 2 ) < 0, 则 x G 2

97 多总体 Fisher 判别 如果变量有多个总体, 通常要选多个投影, 即若干判别函数来进行判别 设有 K 个总体 G 1, G 2,, G k, 它们有共同的协方差阵, 均值向量分别为 μ 1, μ 2,, μ k, 考虑 p 维随机向量 x 的线性组合 y = a T x; y 相对 G i 的均值为 μ iy = E y x G i 方差 :Var y = a T a. 令 μ y = 1 k μ iy = 1 k = a T μ i, i = 1,, k; a T μ i =a T μ, 其中 μ = 1 k μ i 同样考虑比值 i=1 k (μiy μ y ) 2 var(y) = i=1 k (a T μ i a T μ) 2 var(y) = at Ga, 其中 G = a T a i=1 (μ i k μ)(μ i μ) T

98 考察 不妨设 a T a=1 a T Ga a T a 上式变为 b T 1 2G 1 2b, 其中 b = 1 2a,a = 1 2b 也就是当 b 满足单位为 1, 求 b T 1 2G 1 2b 的最大值, 即求在单位球上二次型的最大值 设 λ 1, λ 2,, λ s 为半正定矩阵 1 2G 1 2 的前 s 个不为零的特征根, 并从小到大排序 α 1, α 2,, α s 为对应的单位特征向量,α s+1, α s+2,, α p 对应 0 的单位特征向量 (1) 当上述的 b=α 1 时, 即 a = 1 2α 1, at Ga 达到最值 λ a T a 1, 此时称 y = a T x 为第一判别函数 T (2) 在约束条件 α k+1 α i = 0, i = 1,, k, 称 y = ( 1 2α k+1 ) T x 为第 k+1 个判别函数 S 个判别函数的性质?

99 1 第一判别函数 第二判别函数, 怎样解释? 2 以求极大值和上述方法求矩阵的特征根的区别是什么? 求出的投影方向不一样! R 语言中的函数 library(mass)## 用判别函数 lda() attach(iris)# 把数据变量的名字放入内存, 这样能直接使用各列数据 ld<-lda(species~sepal.length+sepal.width+petal.length+petal.width) pre<-predict(ld) Newclass<-pre$class Cbind(species, newclass,pre$x)

100 Call: lda(species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width) Prior probabilities of groups: setosa versicolor virginica Group means: Sepal.Length Sepal.Width Petal.Length Petal.Width setosa versicolor virginica 输出中包括 lda() 所用公式 先验概率 各组均值向量 第一及第二线性判别函数的系数 两个判别式对区分各总体贡献的大小等 Coefficients of linear discriminants: LD1 LD2 Sepal.Length Sepal.Width Petal.Length Petal.Width Proportion of trace: LD1 LD

101 pre<-predict(ld) Newclass<-pre$class Cbind(Species, newclass,pre$x) 这里 Species,newclass 是回判类别, LD1 和 LD2 分别是第一 第二线性判别函数的值 我们可以用 table() 函数来列表比较 Table(newclass,Species) table(newclass,species)

102 7.3 Bayes 判别法 思想 上述两种判别法的缺点 : 一是判别法与总体各自出现的概率大小无关 ; 二是判别法与错判之后所造成的损失无关 Bayes 判别法就是为了解决这些问题提出的一种判别方法 它假定对研究对象已经有了一定的认识, 这种认识可以用先验概率来描述, 当取得样本后, 就可以利用样本来修正已有的先验概率, 得到后验分布, 再通过后验分布进行统计推断 Bayes 判别法属于概率判别法, 判别准则则是以个体归属某类的概率最大或错判总平均损失最小为标准

103 面积与概率 C B A 考察 p C B p(a B)

104 p(a, B) p A B = p(b) 说明事件 A 在 B 中出现的几率, 也可以说 A 在 B 中表现的特征相似度 p(c, B) p C B = p(b) C 在 B 中出现的几率,C 在 B 中表现的相似度 如果 p C B < p(a B) 判断 B 落在 A 中的概率比在 C 中大, 即判断为 B A.

105 p ω 1 x > p ω 2 x x ω 1 p(x ω 1 ) > p(x ω 2 ) x ω 1

106 我们的手画图表示 p x ω 1, p x ω 2 作为 x 函数的分布图, 中间过相交点的线, 把特征空间 (feature space) 分为两部分区域 R 1 和 R 2. 根据贝叶斯决策理论

107 7.3.1 两总体的 Bayes 判别 设有两个总体 G 1 和 G 2, 它们的概率密度函数分别为 f 1 (x) 与 f 2 x, 其中 x 是一个 p 维随机向量,Ω 为 x 的所有可能取值构成的样本空间,R 1 为 x 的根据某种规则被判入总体 G 1 的取值全体的集合, 那么 R 2 = Ω R 1 就为 x 的根据某种规则被判入总体 G 2 的取值全体的集合 设样本 x 来自总体 G 1 ( 记为 x G 1 ), 但被判为总体 G 2 的概率为 P 2 1 = P x R 2 x G 1 = R2 f 1 x dx 又记 x 来自总体 G 2, 但被判入总体 G 1 的概率 P 1 2 = P x R 1 x G 2 = R1 f 2 x dx

108 类似地, x 来自总体 G 1 被判入 G 1, 来自 G 2 被判入 G 2 概率分别为 : P 1 1 = P x R 1 x G 1 = R 1 f 1 x dx P 2 2 = P x R 2 x G 2 = R2 f 2 x dx. 又设总体 G 1 和 G 2 的先验概率 (prior probabilities) 分别为 p 1 和 p 2, 且 p 1 + p 2 = 1 于是 P x 被正确判入 G 1 = P x G 1, x R 1 = P x R 1 x G 1 P x G 1 = P 1 1. p 1

109 P x 被错误判入 G 1 = P x G 2, x R 1 = P x R 1 x G 2 P x G 2 = P 1 2. p 2 同理 P x 被正确判入 G 2 = P(2 2)p 2 P x 被错误判入 G 2 = P(2 1)p 1 假设 L j i, (i, j = 1,2) 表示 x 来自总体 G i 而被误判入总体 G j 引起的损失, 显然有 L 1 1 = L 2 2 = 0, 将上述误判概率与误判损失结合, 可以定义平均误判损失 (expected cost of misclassification, 简记为 ECM) 为 : ECM R 1, R 2 = L 2 1 P 2 1 p 1 + L 1 2 P(1 2)p 2

110 合理的判别选择为 ECM 极小化 极小化 ECM 对应的样本空间的划分为 R 1 = x f 1 x f 2 x L 1 2 L 2 1. p 2 p 1 R 2 = x f 1 x f 2 x < L 1 2 L 2 1. p 2 p 1 上式称为 Bayes 判别的判别准则

111 命题 : 设 α β 为 p 为列向量, 则 αα T + ββ T 矩阵的秩小于等于 2 证明 : 令 X = αt β T, 则 X T X = αα T + ββ T. 又因为 r X T X = r X, X 为 2 p 矩阵, 所以 r X 2. 则 r X = 2 充分必要条件为 α β 线性无关 同理可以证明 : 设 α i,i = 1,2,, m, 为 m 个 p 为列向量 (m<p), 则 α 1 α 1 T + α 2 α 2 T + + α m α m T 矩阵的秩小于等于 m 秩等于 m 的充分必要条件为它们线性无关

112 对 Fisher 判别法来说, 它的判别函数个数不会高于类的个数 进一步的问题 : 已知几个判别函数, 明确能够分出几类? 判别函数的个数与明确判别的类是否一致? 也就是说判别函数的作用怎样具体体现? 对损失的公式解释, 参见手画图

113 R 函数 :lda() 不同于 Fisher 判别法, 用这个函数时首先输入先验概率, 并假设误判损失相等 library(mass)## 用判别函数 lda() attach(iris)# 把数据变量的名字放入内存, 这样能直接使用各列数据 ld<-lda(species~sepal.length+sepal.width+petal.length+petal.width,, prior=c(1/3,1/3,1/3)) pre<-predict(ld)# 回代判别 newclass<-pre$class cbind(species, newclass,pre$x) table(newclass,species)