R软件应用多元分析演示文稿

R软件应用多元分析演示文稿当前1页，总共52页。R软件应用多元分析当前2页，总共52页。8.1判别法则(分类)已知有多少类，并且在训练样本的前提下，利用训练样本得到判别函数，对待测样本进行分类；8.1.1距离判别判别问题，就是将p维欧几里得空间Rp划分成k个互不相交的区域R1,R2,…,Rk。若x∈Ri,i=1,2,…,k,则判定x属于总体Xi,i=1,2,…,k.Mahalanobis距离的概念：定义8.1设x,y是从均值为μ,协方差矩阵为Σ的总体X中抽取的两个样本,则总体X内两点x,y的Mahalanobis距离定义为样本x与总体X的Mahalanobis距离为：当前3页，总共52页。例如：=1.661.662.34从欧氏距离看A到μ1的距离比到μ2的距离要近，但从概率分布的角度看,,说明A到μ2的距离比到μ1的距离要近.标准化Mahalanobis距离符合概率分布内涵.当前4页，总共52页。2.判别准则与判别函数2.1两个总体的距离判定.总体X1,X2的均值向量分别为μ1,μ2,协方差分别为Σ1,Σ2,给定样本x，判断x来自哪一个总体.1.μ1≠μ2,

Σ1=Σ2判断准则:判断准则:当前5页，总共52页。总体的均值与协方差未知时：设是来自总体X1的n1个样本,是来自总体X2的n2个样本,则样本的均值与协方差阵为判断准则:当前6页，总共52页。1.μ1≠μ2,

Σ1≠Σ2判断函数:总体的均值与协方差未知时：总体的均值与协方差已知时：MahalanobisDistanceReturnsthesquaredMahalanobisdistanceofallrowsinxandthevectormu=centerwithrespecttoSigma=cov.Thisis(forvectorx)definedasD^2=(x-μ)'Σ^-1(x-μ)Usage:mahalanobis(x,center,cov,inverted=FALSE,...)X:vectorormatrixofdatawith,say,pcolumns.

Center:meanvectorofthedistributionorseconddatavectoroflengthp.Cov:covariancematrix(pxp)ofthedistribution.当前7页，总共52页。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))elseif(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))#产生一个行矩阵，共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(iin1:nx){if(w[i]>0)blong[i]<-1elseblong[i]<-2}blong}#X1，X2类的训练样本#TstX=NULL待测样本为2个训练样本之和#数据全部转化成矩阵,行表示样本个数，列表示样本维数n#根据第i个样本的wi值,返回样本类别结果理论中的样本按列排列X=(X1,X2,…,Xn)，每列是一个样本，n列表示n个样本，这里样本按行排X=(X1,X2,…,Xn)T当前8页，总共52页。4.判别实例例8.1在研究砂基液化问题中，选了7个因子，今从已液化和未液化的地层中分别抽了12个和23个样本，数据列在表中，其中I类表示已液化类，II类表示未液化类。试建立距离判别的判别准则，并按判别准则对原35个样本进行回代(即按判别准则进行分类)，分析误判情况。编号类别x1x2x3x4x5x6x71I6.6391660.12202I6.63916120.12203I6.1471660.08124I6.14716120.08125I8.43227.5190.35756I7.2617280.3307I8.41133.56180.15758I7.55216120.16409I7.5523.57.560.164010I8.311307.5350.1218011I7.817213.5140.214512I7.81721.53150.214513II8.4321540.357514II8.43229100.357515II8.4322.54100.357516II6.3114.57.530.21517II784.54.590.253018II7867.540.253019II781.5610.253020II8.31611.5440.087021II8.31610.52.510.087022II7.263.54120.33023II7.261330.33024II7.261650.33025II5.562.5370.181826II8.41133.54.560.157527II8.41133.54.580.157528II7.5521660.164029II7.55217.580.164030II8.3970650.1518031II8.3972.5650.1518032II8.38906100.1618033II8.3561.56130.2518034II7.817213.560.214535II7.828314.560.1845#R里的数据就是这样排,样本均值是对每个指标按列求均值,然后组成样本均值当前9页，总共52页。R实现：classx1=read.table('dataexample801x1.txt')classx2=read.table('dataexample801x2.txt')discriminiant.distance(classx1,classx2,var.equal=T)

123456789101112131415161718192021222324252611111111211122222222222222272829303132333435blong211222222blong#在认为两个总体协方差相同的情况下，有3个点判错discriminiant.distance(classx1,classx2)

123456789101112131415161718192021222324252611111111211122222222222222272829303132333435blong222222222blong#在认为两个总体协方差不同的情况下，有1个点判错当前10页，总共52页。5.多分类问题的距离判别μ1≠μ2…≠μk,

Σ1=Σ2…=Σ

k相应的判别准则：distinguish.distance<-function(TrnX,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))

elseif(is.matrix(TstX)!=TRUE)

#待测样本是多样本,但不是矩阵形式时TstX<-as.matrix(TstX)#转成矩阵(如data.frame类型转成矩阵)if(is.matrix(TrnX)!=TRUE)TrnX<-as.matrix(TrnX)nx<-nrow(TstX)blong<-matrix(rep(0,nx),nrow=1,dimnames=list(“blong”,1:nx))#本页语句都是准备工作#如果TrnG从主函数未接收到因子数据#待测样本TstX是单样本时候,是向量vector,此时将其转为矩阵(是列矩阵),然后再转成行矩阵#则是2分类问题,而非多分类,可省略#行名称为”blong”,列名称为数字1到nx#产生类别矩阵blong,初始值全为0当前11页，总共52页。Continue:g<-length(levels(TrnG))mu<-matrix(0,nrow=g,ncol=ncol(TrnX))

for(iin1:g)mu[i,]<-colMeans(TrnX[TrnG==i,])

D<-matrix(0,nrow=g,ncol=nx)#ncol个样本因子按列排,g个类别按行排#对属于第i个类的样本求他们因子的均值,结果存到mu的第i行#产生0阵,行数为类别数g,列数为样本数nx#得到多分类的类别,共g个:[1]2[,1][,2][,3][,4][,5][,6][,7][1,]0000000[2,]0000000[,1][,2][,3][,4][,5][,6][,7][1,]7.35833373.666671.4583336.0000015.2500000.171666749.50000[2,]7.68695769.608702.0434785.239136.3478260.215652270.34783

[,1][,2][,3][,4][,5][,6][,7][,8][,9][,10][,11][,12][,13][1,]0000000000000[2,]0000000000000[,26][,27][,28][,29][,30][,31][,32][,33][,34][,35][1,]0000000000[2,]0000000000当前12页，总共52页。Continue:if(var.equal==TRUE||var.equal==T){for(iin1:g)D[i,]<-mahalanobis(TstX,mu[i,],var(TrnX))}else{for(iin1:g)D[i,]<-mahalanobis(TstX,mu[i,],var(TrnX[TrnG==i,]))}for(jin1:nx){dmin<-Inffor(iin1:g)if(D[i,j]<dmin){dmin<-D[i,j];blong[j]<-i}}blong#待测样本到第i类的马氏距离

[,1][,2][,3][,4][,5][,6][,7][,8][1,]181.3889182.5306162.9359164.1592233.7561205.8525238.8812214.1178[2,]181.3889182.5306162.9359164.1592233.7561205.8525238.8812214.1178[,9][,10][,11][,12][,13][,14][,15][,16][1,]219.4913201.6875185.8174184.6754222.2703241.7303218.1065169.8246[2,]219.4913201.6875185.8174184.6754222.2703241.7303218.1065169.8246[,17][,18][,19][,20][,21][,22][,23][,24][1,]182.1071197.3248186.0215224.1768221.5524185.4759183.1583192.0052[2,]182.1071197.3248186.0215224.1768221.5524185.4759183.1583192.0052[,25][,26][,27][,28][,29][,30][,31][,32][1,]114.8461229.1703229.5824213.1659220.9516181.5381180.7470181.3109[2,]114.8461229.1703229.5824213.1659220.9516181.5381180.7470181.3109[,33][,34][,35][1,]175.8290184.0261181.7853[2,]175.8290184.0261181.7853#对第j个样本,纵向求min，如果该最小值位于第i行，则第j个样本就是属于第i类方差未知方差已知当前13页，总共52页。8.1.2Bayes判别1.误判概率与误判损失x被判为X2x实际来自X1来自X2,但被判为x1的概率:来自X1,但被判为x2的概率:来自X1,但被判为x1的概率:来自X2,但被判为x2的概率:总体X1的先验概率平均误判损失ECM:ECM(R1,R2)=L(2|1)P(2|1)p1+L(1|2)P(1|2)p2来自X1被判为X2引起的损失来自X2被判为X1引起的损失当前14页，总共52页。2.两个总体的Bayes判别ECM(R1,R2)=L(2|1)P(2|1)p1+L(1|2)P(1|2)p20ECM=min划分区域R1和R2:作为Bayes判别准则须计算当前15页，总共52页。正态分布的情况:Xi~N(μi,∑i)(i=1,2)1.∑1=∑2类似地,2.∑1≠∑2当前16页，总共52页。3.R程序与例子R程序略；例8.3下表是某气象站预报有无春旱的实际资料，x1与x2是综合预报因子，有春旱的是6个年份的资料，无春旱的是8个年份的资料，它们的先验概率分别用6/14和8/14来估计，并假设误判损失相等，试用Bayes估计对数据进行分析。序号春旱无春旱124.8-222.1-0.7224.1-2.421.6-1.4326.6-322-0.8423.5-1.922.8-1.6525.5-2.122.7-1.5627.4-3.121.5-1722.1-1.2821.4-1.3当前17页，总共52页。R实现x1=scan('dataexample803x1.txt')x2=scan('dataexample803x2.txt')dim(x1)=c(2,6)[,1][,2][,3][,4][,5][,6][1,]24.824.126.623.525.527.4[2,]-2.0-2.4-3.0-1.9-2.1-3.1dim(x2)=c(2,8)x1=t(x1)x2=t(x2)source('discriminiant.bayes.R')discriminiant.bayes(x1,x2,rate=8/6,var.equal=T)

[,1][,2][1,]24.823.5[2,]-2.0-1.9[3,]24.125.5[4,]-2.4-2.1[5,]26.627.4[6,]-3.0-3.1

1234567891011121314blong11121122222222#4号样本被错判discriminiant.bayes(x1,x2,rate=8/6)

1234567891011121314blong11111122222222#无错判当前18页，总共52页。４.多分类问题的Bayes判别样本共分k类:X1,X2,…,Xk,相应的先验概率为p1,p2,…,pk,假定所有的错判损失相同,则判别准则为:1.∑1=…=∑k=∑2.∑1≠…≠

∑k当前19页，总共52页。R程序……if(var.equal==TRUE||var.equal==T){for(iin1:g){d2<-mahalanobis(TstX,mu[i,],var(TrnX))D[i,]<-d2-2*log(p[i])}}else{for(iin1:g){S<-var(TrnX[TrnG==i,])d2<-mahalanobis(TstX,mu[i,],S)D[i,]<-d2-2*log(p[i])-log(det(S))}}for(jin1:nx){dmin<-Inffor(iin1:g)if(D[i,j]<dmin){dmin<-D[i,j];blong[j]<-i}}blong}当前20页，总共52页。例8.4用Bayes判别对FisherIris数据进行分析.假设先验概率相同,均为1.考虑总体协方差阵不同的情况.x=iris[,1:4]g=gl(3,50)distinguish.bayes(x,g)

123456789101112131415161718192021222324252627blong111111111111111111111111111282930313233343536373839404142434445464748495051blong111111111111111111111112525354555657585960616263646566676869707172737475blong222222222222222223232322767778798081828384858687888990919293949596979899blong223222223222222222222222100101102103104105106107108109110111112113114115116117blong233333333333333333118119120121122123124125126127128129130131132133134135blong333333333333333333136137138139140141142143144145146147148149150blong333333333333333#误判概率为1-145/150=3.33%当前21页，总共52页。8.1.3Fisher判别按类内方差尽量小,类间方差尽量大的准则求判别函数.(以2个总体为例)判别准则总体X1,X2的均值与协方差阵分别为μ1,μ2和Σ1,Σ2,对于样本x,考虑其判别函数:判别准则为:U(x)=?当前22页，总共52页。2.线性判别函数中系数的确定u(x)为线性函数设总体X1,X2的样本容量为n1,n2;则u1,u2和σ1,σ

2的估计:当前23页，总共52页。确定判别函数

若：进一步：判别准则为:当前24页，总共52页。4.R程序与例子例8.5用Fisher判别解例8.1classx1=read.table('dataexample801x1.txt')classx2=read.table('dataexample801x2.txt')discriminiant.fisher(classx1,classx2)结果：

1234567891011121314151617181920212223242526blong111111111111222222222222222728

29303132333435blong211222222#28，29号样本为误判样本R程序…

mu1<-colMeans(TrnX1);mu2<-colMeans(TrnX2)

S<-(n1-1)*var(TrnX1)+(n2-1)*var(TrnX2)

mu<-n1/(n1+n2)*mu1+n2/(n1+n2)*mu2w<-(TstX-rep(1,nx)%o%mu)%*%solve(S,mu2-mu1);…当前25页，总共52页。8.2聚类分析[7]常用的几种距离：第2个样本与第n个样本之间的距离记为d2nordn28.2.1距离和相似系数绝对值距离or“棋盘距离”or“城市街区”距离Euclid(欧几里得)距离Minkowski(闵可夫斯基)距离当前26页，总共52页。continueChebyshev(切比雪夫)距离Mahalanobis距离为：Lance和Williams距离定性变量样本间的距离第i个样本记为:项目项目的类目数样本x(1)x(2)性别

外语 专业 职业男 女英 日德俄 统计会计金融教师 工程师1010000010101100010010类目当前27页，总共52页。Continue:样本x(1)x(2)性别 外语 专业 职业男 女英 日德俄 统计会计金融教师 工程师10100000101011000100101-1配对0-0配对不配对第i个样本和第j个样本在第k个项目的第l类上1-1配对第i个样本和第j个样本在第k个项目的第l类上0-0配对第i个样本和第j个样本在第k个项目的第l类上不配对0-0配对数1-1配对数不配对数表中的样本距离d12=6/7=0.8571429当前28页，总共52页。R中的距离函数Usagedist(x,method="euclidean",diag=FALSE,upper=FALSE,p=2)DescriptionThisfunctioncomputesandreturnsthedistancematrixcomputedbyusingthespecifieddistancemeasuretocomputethedistancesbetweentherowsofadatamatrix.x:anumericmatrix,dataframeor"dist"object.method:thedistancemeasuretobeused.Thismustbeoneof"euclidean","maximum","manhattan","canberra","binary"or"minkowski".Anyunambiguoussubstringcanbegiven.Euclidean:euclid距离Maximum:chebyshev距离Manhattan:绝对值距离Canberra:lance距离Binary:定性变量的距离Minkowski:minkowski距离当前29页，总共52页。2.数据中心化与标准化表示因素表示样本对第i个样本中心化，标准化当前30页，总共52页。(1)极差标准化.对第i个样本极差标准化,正规化x=data.frame(x1=c(1,2),x2=c(1,2),x3=c(1,2),x4=c(1,2))x1x2x3x41111122222ap=apply(x,2,mean)x1x2x3x41.51.51.51.5center=sweep(x,2,ap)r=apply(x,2,max)-apply(x,2,min)x_star=sweep(center,2,r,'/')x_star=sweep(center,2,sd(x),'/')>rx1x2x3x41111>centerx1x2x3x41-0.5-0.5-0.5-0.520.50.50.50.5>x_starx1x2x3x41-0.5-0.5-0.5-0.520.50.50.50.5>x_starx1x2x3x41-0.7071068-0.7071068-0.7071068-0.707106820.70710680.70710680.70710680.7071068#普通标准化当前31页，总共52页。continue(2)极差正规化变换.x=data.frame(x1=c(1,2),x2=c(1,2),x3=c(1,2),x4=c(1,2))x1x2x3x41111122222ap=apply(x,2,min)x1x2x3x41111center=sweep(x,2,ap)r=apply(x,2,max)-apply(x,2,min)x_star=sweep(center,2,r,'/')>rx1x2x3x41111>centerx1x2x3x41000021111>x_starx1x2x3x41000021111当前32页，总共52页。3.相似系数相似系数用于对变量进行分类。夹角余弦Xi与Xj的夹角余弦称为两向量的相似系数x=data.frame(x1=c(1,2),x2=c(1,2),x3=c(1,2),x4=c(1,2))y=scale(x,center=F,scale=T)/sqrt(nrow(x)-1)x1x2x3x4[1,]0.44721360.44721360.44721360.4472136[2,]0.89442720.89442720.89442720.8944272c=t(y)%*%y将样本列(变量)标准化当前33页，总共52页。相关系数中心化样本(变量)的相关矩阵R实现：cor(x)当前34页，总共52页。8.2.2系统聚类法记号：dij:第i个样本与第j个样本的距离.G1,G2:表示类.DKL:GK与GL的(类)距离.最短距离法类与类之间的距离为两类最近样本间的距离:当某步骤类GK和GL合并为GM后,按最短距离法计算新类GM与其他类GJ的类间距离:最长距离法递推公式当前35页，总共52页。中间距离法推广：类平均法可变类平均法:类GK的样本个数Mcquitty相似分析递推公式当前36页，总共52页。类与类之间的距离定义为他们重心(均值)之间的Euclid距离.设GK和GL的重心分别为和.重心法递推公式离差平方和法(ward方法)递推公式GK和GL的平方距离也可定义为：与重心法相差一个系数，表明表明大样本类不易合并，这更符合实际。当前37页，总共52页。7.R相关函数及其用法Usagehclust(d,method="complete",members=NULL)DescriptionHierarchicalclusteranalysisonasetofdissimilaritiesandmethodsforanalyzingit.d:adissimilaritystructureasproducedbydist.method:"ward","single","complete","average","mcquitty","median"or"centroid".例8.6设有5个样本,每个样本只有一个指标,分别是1,2,6,8,11,样本间的距离选用Euclid距离,试用最短距离法、最长距离法等方法进行聚类分析，并画出相应的谱系图.R实现:x=c(1,2,6,8,11)dim(x)=c(5,1)d=dist(x)>d

1234213544762510953第一个样本到第2,3,4,5个样本的距离当前38页，总共52页。例8.6(续)hc1=hclust(d,'single')hc2=hclust(d,'complete')hc3=hclust(d,'median')hc4=hclust(d,'mcquitty')hc5=hclust(d,'average')hc6=hclust(d,'centroid')hc7=hclust(d,'ward')opar=par(mfrow=c(3,3))plot(hc1,hang=-1)plot(hc2,hang=-1)plot(hc3,hang=-1)plot(hc4,hang=-1)plot(hc5,hang=-1)plot(hc6,hang=-1)plot(hc7,hang=-1)当前39页，总共52页。例8.7对305名女中学生测量8个体型指标,相应的相关矩阵如表，将相关系数看成相似系数,定义距离为dij=1-rij,用最长距离法做系统分析.身高手臂长上肢长下肢长体重颈围胸围胸宽x1x2x3x4x5x6x7x8x11x20.8461x30.8050.8811x40.8590.8260.8011x50.4730.3760.380.4361x60.3980.3260.3190.3290.7621x70.3010.2770.2370.3270.730.5831x80.3820.2770.3450.3650.6290.5770.5391当前40页，总共52页。R实现x=scan('dataexample807.txt')r=as.matrix(x)dim(r)=c(8,8)d=as.dist(1-r)hc=hclust(d)plot(hc,hang=-1)>d123456720.15430.1950.11940.1410.1740.19950.5270.6240.6200.56460.6020.6740.6810.6710.23870.6990.7230.7630.6730.2700.41780.6180.7230.6550.6350.3710.4230.461当前41页，总共52页。8.类个数的确定给定一个阈值.观测样本的散点图.(仅限于二维,三维样本)试用统计量.根据谱系图确定分类个数的原则:A.各类重心的距离必须很大.B.确定的类中,各类所包含的元素都不要太多.C.类的个数必须符合实用的目的.D.若采用不同的聚类方法处理,则在各自的聚类图中应发现相同的类.Usagerect.hclust(tree,k=NULL,which=NULL,x=NULL,h=NULL,border=2,cluster=NULL)DescriptionDrawsrectanglesaroundthebranchesofadendrogramhighlightingthecorrespondingclusters.Tree是由hclust生成的对象;K是类的个数;H是谱系图中的阈值;当前42页，总共52页。Rect.hclust()实例(8.7)当前43页，总共52页。9.实例表中给出了1999年全国31个省,市,自治区的城镇居民家庭平均每人全年消费性支出的8个主要指标(变量)数据.这8个变量是:x1:食品;x2:衣着;x3:家庭设备用品及服务;x4:医疗保障;x5:交通与通信;x6:娱乐教育文化服务;x7:居住;x8:杂项商品和服务;分别使用最长距离法,类平均法,重心法和ward方法对各地区做聚类分析.x1x2x3x4x5x6x7x8北京2959.19730.79749.41513.34467.871141.82478.42457.64天津2459.77495.47697.33302.87284.19735.97570.84305.08河北1495.63515.9362.37285.32272.95540.58364.91188.63山西1046.33477.77290.15208.57201.5414.72281.84212.1内蒙1303.97524.29254.83192.17249.81463.09287.87192.96辽宁1730.84553.9246.91279.81239.18445.2330.24163.86吉林1561.86492.42200.49218.36220.69459.62360.48147.76黑龙江1410.11510.71211.88277.11224.65376.82317.61152.85上海3712.31550.74893.37346.935271034.98720.33462.03浙江2629.16557.32689.73435.69514.66795.87575.76323.36安徽1844.78430.29271.28126.33250.56513.18314151.39福建2709.46428.11334.12160.77405.14461.67525.13232.29江西1563.78303.65233.81107.9209.7393.99509.39160.12当前44页，总共52页。continue山东1675.75613.32550.71219.79272.59599.43371.62211.84河南1427.65431.79288.55208.14217337.76421.31165.32湖北1783.43511.88282.84201.01237.6617.74523.52182.52湖南1942.23512.27401.39206.06321.29697.22492.6226.45广东3055.17353.23564.56356.27811.88873.061082.82420.81广西2033.87300.82338.65157.78329.06621.74587.02218.27海南2057.86186.44