R软件应用多元分析

R软件应用多元分析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距离为：例如:=1、661.662.34从欧氏距离看A到μ1得距离比到μ2得距离要近,但从概率分布得角度看,,说明A到μ2得距离比到μ1得距离要近、标准化Mahalanobis距离符合概率分布内涵.2、判别准则与判别函数2、1两个总体得距离判定、总体X1,X2得均值向量分别为μ1,μ2,协方差分别为Σ1,Σ2,给定样本x,判断x来自哪一个总体、1、μ1≠μ2,

Σ1=Σ2判断准则:判断准则:总体得均值与协方差未知时:设是来自总体X1的n1个样本,是来自总体X2的n2个样本,则样本的均值与协方差阵为判断准则: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、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)T4、判别实例例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里的数据就是这样排,样本均值是对每个指标按列求均值,然后组成样本均值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个点判错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,初始值全为0Continue: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,]0000000000Continue: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类方差未知方差已知大家有疑问的，可以询问和交流可以互相讨论下，但要小声点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引起的损失2、两个总体得Bayes判别ECM(R1,R2)=L(2|1)P(2|1)p1+L(1|2)P(1|2)p20ECM=min划分区域R1和R2:作为Bayes判别准则须计算正态分布得情况:Xi~N(μi,∑i)(i=1,2)1、∑1=∑2类似地,2、∑1≠∑23、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、3R实现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#无错判４、多分类问题得Bayes判别样本共分k类:X1,X2,…,Xk,相应得先验概率为p1,p2,…,pk,假定所有得错判损失相同,则判别准则为:1、∑1=…=∑k=∑2、∑1≠…≠

∑kR程序……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}例8、4用Bayes判别对FisherIris数据进行分析、假设先验概率相同,均为1、考虑总体协方差阵不同得情况、x=iris[,1:4]g=gl(3,50)distinguish、bayes(x,g)

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

2得估计:确定判别函数

若：进一步:判别准则为:4、R程序与例子例8、5用Fisher判别解例8、1classx1=read、table('dataexample801x1、txt')classx2=read、table('dataexample801x2、txt')discriminiant、fisher(classx1,classx2)结果:

1234567891011121314151617181920212223242526blong11111111111122222222222222272829303132333435blong211222222#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);…8、2聚类分析[7]常用得几种距离:第2个样本与第n个样本之间的距离记为d2nordn28、2、1距离和相似系数绝对值距离or“棋盘距离”or“城市街区”距离Euclid(欧几里得)距离Minkowski(闵可夫斯基)距离continueChebyshev(切比雪夫)距离Mahalanobis距离为：Lance和Williams距离定性变量样本间得距离第i个样本记为:项目项目的类目数样本x(1)x(2)性别 外语 专业 职业男 女英 日德俄 统计会计金融教师 工程师1010000010101100010010类目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、8571429R中得距离函数Usagedist(x,method="euclidean",diag=FALSE,upper=FALSE,p=2)DescriptionThisfunctionputesandreturnsthedistancematrixputedbyusingthespecifieddistancemeasuretoputethedistancesbetweentherowsofadatamatrix、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距离2、数据中心化与标准化表示因素表示样本对第i个样本中心化,标准化(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#普通标准化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_starx1x2x3x410000211113、相似系数相似系数用于对变量进行分类。夹角余弦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将样本列(变量)标准化相关系数中心化样本(变量)得相关矩阵R实现:cor(x)8、2、2系统聚类法记号:dij:第i个样本与第j个样本得距离、G1,G2:表示类、DKL:GK与GL得(类)距离、最短距离法类与类之间的距离为两类最近样本间的距离:当某步骤类GK和GL合并为GM后,按最短距离法计算新类GM与其她类GJ得类间距离:最长距离法递推公式中间距离法推广：类平均法可变类平均法:类GK的样本个数Mcquitty相似分析递推公式类与类之间的距离定义为他们重心(均值)之间的Euclid距离.设GK和GL的重心分别为和.重心法递推公式离差平方和法(ward方法)递推公式GK和GL的平方距离也可定义为：与重心法相差一个系数,表明表明大样本类不易合并,这更符合实际。7、R相关函数及其用法Usagehclust(d,method="plete",members=NULL)DescriptionHierarchicalclusteranalysisonasetofdissimilaritiesandmethodsforanalyzingit、d:adissimilaritystructureasproducedbydist、method:"ward","single","plete","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)>d1234213544762510953第一个样本到第2,3,4,5个样本的距离例8、6(续)hc1=hclust(d,'single')hc2=hclust(d,'plete')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)例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、5391R实现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、4618、类个数得确定给定一个阈值、观测样本得散点图、(仅限于二维,三维样本)试用统计量、根据谱系图确定分类个数得原则:A、各类重心得距离必须很大、B、确定得类中,各类所包含得元素都不要太多、C、类得个数必须符合实用得目得、D、若采用不同得聚类方法处理,则在各自得聚类图中应发现相同得类、Usagerect、hclust(tree,k=NULL,which=NULL,x=NULL,h=NULL,border=2,cluster=NULL)DescriptionDrawsrectanglesaroundthebranchesofadendrogramhighlightingthecorrespondingclusters、Tree就是由hclust生成得对象;K就是类得个数;H就是谱系图中得阈值;Rect、hclust()实例(8、7)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、12continue山东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、56