田茂再《多元统计分析 》书中代码汇总

《多元统计分析》代码第1章多元分布1.5样本分布和极限定理例1.8n=5#n=50X=matrix(0,n,1000)for(iin1:1000){X[,i]=rbinom(n,1,1/2)}XX=(colMeans(X)-1/2)*5^(1/2)plot(density(XX),type="l",xlab="1000Random.Sample",ylab="EstimatedandnormalDensity",main="AsymptoticDistribution,N=5",col="blue",xlim=c(-2,2))par(new=TRUE)dat=rnorm(100000,0,1/2)plot(density(dat),col="red",main="",axes=FALSE)1.6厚尾分布curve(dnorm(x),-5,5,xlab="X",ylab="Y",main="distributioncomparison",xaxp=c(-5,5,2),ylim=c(0,0.4),col="blue")curve(dcauchy(x,0,1),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.4),col="red")legend("topright",c("Gaussian","Cauchy"),pch=c(16,16),col=c("blue","red"),text.col=c("blue","red"),bty="n")abline(v=c(-2,-1,1,2))text(-2.1,0.39,"-2f",col="blue",adj=c(0,-0.1))text(-1.1,0.39,"-f",col="blue",adj=c(0,-0.1))text(0.9,0.39,"f",col="blue",adj=c(0,-0.1))text(1.9,0.39,"2f",col="blue",adj=c(0,-0.1))1.6.1广义双曲分布library(Runuran)#为广义双曲分布产生分布函数目标distr1<-udghyp(lambda=0.5,alpha=1,beta=0,delta=1,mu=0)#制造生成器目标；利用PINV(逆)方法gen1<-pinvd.new(distr1)#产生大小为10000的样本x1<-ur(gen1,10000)#概率密度函数dx1<-ud(gen1,x1)#累积分布函数fx1<-up(gen1,x1)o1=order(x1)#为HYD产生分布目标distr2<-udghyp(lambda=1,alpha=1,beta=0,delta=1,mu=0)gen2<-pinvd.new(distr2)x2<-ur(gen2,10000)dx2<-ud(gen2,x2)fx2<-up(gen2,x2)o2=order(x2)#画图op<-par(mfrow=c(1,2))plot(sort(x1),dx1[o1],type="l",xlab="X",ylab="Y",main="pdfofGH,HYPandNIC",xlim=c(-5,5),xaxp=c(-5,5,2),ylim=c(0,0.52),col="black")lines(sort(x2),dx2[o2],type="l",col="red")lines(sort(x3),dx2[o3],type="l",col="blue")legend("topright",c("GH","NIC","HYP"),pch=c(16,16,16),col=c("black","red","blue"),text.col=c("black","red","blue"),bty="n")plot(sort(x1),fx1[o1],type="l",xlab="X",ylab="Y",main="cdfofGH,HYPandNIC",xlim=c(-5,5),xaxp=c(-5,5,2),ylim=c(0,1),col="black")lines(sort(x2),fx2[o2],type="l",col="red")lines(sort(x3),fx2[o3],type="l",col="blue")legend("topleft",c("GH","NIG","HYP"),pch=c(16,16,16),col=c("black","red","blue"),text.col=c("black","red","blue"),bty="n")par(op)1.6.2学生t分布图1-5op<-par(mfrow=c(1,2))curve(dt(x,3),-5,5,xlab="X",ylab="Y",main="pdfoft-distribution",xaxp=c(-5,5,2),ylim=c(0,0.42),col="black")curve(dt(x,6),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.42),col="blue")curve(dt(x,30),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.42),col="red")legend("topright",c("t3","t6","t30"),pch=c(16,16,16),col=c("black","blue","red"),text.col=c("black","blue","red"),bty="n")curve(pt(x,3),-5,5,xlab="X",ylab="Y",main="cdfoft-distribution",xaxp=c(-5,5,2),ylim=c(0,1),col="black")curve(pt(x,6),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,1),col="blue")curve(pt(x,30),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,1),col="red")legend("topleft",c("t3","t6","t30"),pch=c(16,16,16),col=c("black","blue","red"),text.col=c("black","blue","red"),bty="n")par(op)图1-6op<-par(mfrow=c(1,2))curve(dt(x,1),2.6,3.8,xlab="X",ylab="Y",main="tailcomparison-t-distribution",xaxp=c(3,3.5,1),ylim=c(0,0.04),col="black")curve(dt(x,3),2.6,3.8,add=TRUE,xaxp=c(3,3.5,1),ylim=c(0,0.04),col="blue")curve(dt(x,9),2.6,3.8,add=TRUE,xaxp=c(3,3.5,1),ylim=c(0,0.04),col="red")curve(dt(x,45),2.6,3.8,add=TRUE,xaxp=c(3,3.5,1),ylim=c(0,0.04),col="deeppink")curve(dnorm(x),2.6,3.8,add=TRUE,xaxp=c(3,3.5,1),ylim=c(0,0.04),col="grey")points(x=3.2,y=dt(3.2,1),cex=0.8,pch=1,col="black")text(x=3.21,y=dt(3.2,1),"t1",cex=0.8,col="black",adj=c(0,-0.1))points(x=3.2,y=dt(3.2,3),cex=1,pch=1,col="blue")text(x=3.21,y=dt(3.2,3),"t3",cex=0.8,col="blue",adj=c(0,-0.1))points(x=3.2,y=dt(3.2,9),cex=1,pch=1,col="red")text(x=3.21,y=dt(3.2,9),"t9",cex=0.8,col="red",adj=c(0,-0.1))points(x=3.2,y=dt(3.2,45),cex=1,pch=1,col="deeppink")text(x=3.21,y=dt(3.2,45),"t45",cex=0.8,col="deeppink",adj=c(0,-0.1))points(x=3.2,y=dnorm(3.2),cex=1,pch=1,col="grey")text(x=3.21,y=0,"Gaussian",cex=0.8,col="grey",adj=c(0,-0.1))curve((abs(x))^(-2),1.2,1.5,xlab="X",ylab="Y",main="tailcomparison-approximation",xaxp=c(1.2,1.5,6),ylim=c(0,0.7),col="black")curve((abs(x))^(-4),1.2,1.5,add=TRUE,xaxp=c(1.2,1.5,6),ylim=c(0,0.7),col="black")curve((abs(x))^(-10),1.2,1.5,add=TRUE,xaxp=c(1.2,1.5,6),ylim=c(0,0.7),col="red")points(x=1.35,y=(abs(1.36))^(-2),cex=0.8,pch=1,col="black")text(x=1.36,y=(abs(1.35))^(-2),"t1",cex=0.8,col="black",adj=c(0,-0.1))points(x=1.35,y=(abs(1.36))^(-4),cex=1,pch=1,col="blue")text(x=1.36,y=(abs(1.35))^(-2),"t3",cex=0.8,col="blue",adj=c(0,-0.1))points(x=1.35,y=(abs(1.36))^(-10),cex=1,pch=1,col="red")text(x=1.36,y=(abs(1.35))^(-10),"t9",cex=0.8,col="red",adj=c(0,-0.1))par(op)1.6.3拉普拉斯分布library(Runuran)#为拉普拉斯分布产生分布函数目标distr1<-udlaplace(location=0,scale=1)#制造生成器目标；利用PINV(逆)方法gen1<-pinvd.new(distr1)#产生大小为10000的样本x1<-ur(gen1,10000)#密度dx1<-ud(gen1,x1)#累积分布函数fx1<-up(gen1,x1)o1=order(x1)distr2<-udlaplace(location=0,scale=1.5)gen2<-pinvd.new(distr2)x2<-ur(gen2,10000)dx2<-ud(gen2,x2)fx2<-up(gen2,x2)o2=order(x2)distr3<-udlaplace(location=0,scale=2)gen3<-pinvd.new(distr3)x3<-ur(gen3,10000)dx3<-ud(gen3,x3)fx3<-up(gen3,x3)o3=order(x3)#绘图op<-par(mfrow=c(1,2))plot(sort(x1),dx1[o1],type="l",xlab="X",ylab="Y",main="pdfofLaplacedistriction",xlim=c(-5,5),xaxp=c(-5,5,2),ylim=c(0,0.52),col="black")lines(sort(x2),dx2[o2],type="l",col="blue")lines(sort(x3),dx2[o3],type="l",col="red")legend("topright",c("L1","L1.5","L2"),pch=c(16,16,16),col=c("black","blue","red"),text.col=c("black","blue","red"),bty="n")plot(sort(x1),fx1[o1],type="l",xlab="X",ylab="Y",main="cdfofLaplacedistribution",xlim=c(-5,6),xaxp=c(-5,5,2),ylim=c(0,1),col="black",pch=20)lines(sort(x2),fx2[o2],type="l",col="blue")lines(sort(x3),fx2[o3],type="l",col="red")legend("topleft",c("L1","L1.5","L2"),pch=c(16,16,16),col=c("black","blue","red"),text.col=c("black","blue","red"),bty="n")par(op)1.6.4柯西分布op<-par(mfrow=c(1,2))curve(dcauchy(x,0,1),-5,5,xlab="X",ylab="Y",main="pdfofCauchydistribution",xaxp=c(-5,5,2),ylim=c(0,0.32),col="black")curve(dcauchy(x,0,1.5),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.32),col="blue")curve(dcauchy(x,0,2),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.32),col="red")legend("topright",c("C1","C1.5","C2"),pch=c(16,16,16),col=c("black","blue","red"),text.col=c("black","blue","red"),bty="n")curve(pcauchy(x,0,1),-5,5,xlab="X",ylab="Y",main="cdfofCauchydistribution",xaxp=c(-5,5,2),ylim=c(0.09,0.91),col="black")curve(pcauchy(x,0,1.5),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0.09,0.91),col="blue")curve(dcauchy(x,0,2),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0.09,0.91),col="red")legend("topleft",c("C1","C1.5","C2"),pch=c(16,16,16),col=c("black","blue","red"),text.col=c("black","blue","red"),bty="n")par(op)1.65混合模型op<-par(mfrow=c(1,2))curve(0.8*dnorm(x)+0.2*dnorm(x,0,3),-9,9,xlabb="X",ylab="Y",main="pdfofGaussianmixtureandGaussian",xaxp=c(-5,5,2),ylim=c(0,0.35),yaxp=c(0,0.3,3),col="blue")curve(dnorm(x,0,1.5),-9,9,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.35),yaxp=c(0,0.3,3),col="red")legend("topright",c("Gaussianmixture","Gaussian"),pch=c(16,16),col=c("blue","red"),text.col=c("blue","red"),bty="n")curve(0.8*pnorm(x)+0.2*pnorm(x,0,3),-9,9,xlabb="X",ylab="Y",main="cdfofGaussianmixtureandGaussian",xaxp=c(-5,5,2),ylim=c(0,1),yaxp=c(0,1,2),col="blue")curve(pnorm(x,0,1.5),-9,9,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,1),yaxp=c(0,1,2),col="red")legend("topleft",c("Gaussianmixture","Gaussian"),pch=c(16,16),col=c("blue","red"),text.col=c("blue","red"),bty="n")par(op)1.6.10广义双曲分布图1-11library(Runuran)#为广义双曲(GH,lambda=1.5)产生分布目标distr1<-udghyp(lambda=1.5,alpha=1,beta=0,delta=1,mu=0)#制造生成器目标；利用PINV(逆)方法gen1<-pinvd.new(distr1)#产生大小为10000的样本x1<-ur(gen1,10000)#概率密度函数dx1<-ud(gen1,x1)#累积分布函数fx1<-up(gen1,x1)o1=order(x1)#为HYD产生分布目标distr2<-udghyp(lambda=1,alpha=1,beta=0,delta=1,mu=0)gen2<-pinvd.new(distr2)x2<-ur(gen2,10000)dx2<-ud(gen2,x2)fx2<-up(gen2,x2)o2=order(x2)#为GH(lambda=0.5)产生分布目标distr3<-udghyp(lambda=0.5,alpha=1,beta=0,delta=1,mu=0)gen3<-pinvd.new(distr3)x3<-ur(gen3,10000)dx3<-ud(gen3,x3)fx3<-up(gen3,x3)o3=order(x3)#为NIG产生分布目标distr4<-udghyp(lambda=-0.5,alpha=1,beta=0,delta=1,mu=0)gen4<-pinvd.new(distr4)x4<-ur(gen4,10000)dx4<-ud(gen4,x4)fx4<-up(gen4,x4)o4=order(x4)#画图op<-par(mfrow=c(1,2))plot(sort(x1),dx1[o1],type="l",xlab="X",ylab="Y",main="tailcomparison-GH",xlim=c(4,5),xaxp=c(4,5,2),ylim=c(0,0.02),col="black")lines(sort(x2),dx2[o2],type="l",col="red")lines(sort(x3),dx2[o3],type="l",col="blue")lines(sort(x4),dx2[o4],type="l",col="deeppink")points(x=4.5,y=ud(gen1,4.5),cex=0.8,pch=1,col="black")text(x=4.55,y=ud(gen1,4.55),"GH(lambda=1.5)",cex=0.8,col="black",adj=c(0,-0.1))points(x=4.5,y=ud(gen2,4.5),cex=1,pch=1,col="blue")text(x=4.55,y=ud(gen2,4.55),"HYP",cex=0.8,col="blue",adj=c(0,-0.1))points(x=4.5,y=ud(gen3,4.5),cex=1,pch=1,col="red")text(x=4.55,y=ud(gen3,4.55),"GH(lambda=0.5)",cex=0.8,col="red",adj=c(0,-0.1))points(x=4.5,y=ud(gen4,4.5),cex=1,pch=1,col="deeppink")text(x=4.55,y=ud(gen4,4.55),"NIG",cex=0.8,col="deeppink",adj=c(0,-0.1))plot(sort(x1),(x1^0.5*exp(-x1))[o1],type="l",xlab="X",ylab="Y",main="tailcomparision-approximation",xlim=c(4,5),xaxp=c(4,5,2),ylim=c(0,0.04),col="black")lines(sort(x2),(exp(-x2))[o2],type="l",col="red")lines(sort(x3),(x3^(-0.5)*exp(-x3))[o3],type="l",col="blue")lines(sort(x4),(x4^(-1.5)*exp(-x4))[o4],type="l",col="deeppink")points(x=4.5,y=4.5^0.5*exp(-4.5),cex=0.8,pch=1,col="black")text(x=4.55,y=4.55^0.5*exp(-4.55),"GH(lambda=1.5)",cex=0.8,col="black",adj=c(0,-0.1))points(x=4.5,y=exp(-4.5),cex=1,pch=1,col="blue")text(x=4.55,y=exp(-4.55),"HYP",cex=0.8,col="blue",adj=c(0,-0.1))points(x=4.5,y=4.5^(-0.5)*exp(-4.5),cex=1,pch=1,col="red")text(x=4.55,y=4.55^(-0.5)*exp(-4.55),"GH(lambda=0.5)",cex=0.8,col="red",adj=c(0,-0.1))points(x=4.5,y=4.5^(-1.5)*exp(-4.5),cex=1,pch=1,col="deeppink")text(x=4.55,y=4.55^(-1.5)*exp(-4.55),"NIG",cex=0.8,col="deeppink",adj=c(0,-0.1))par(op)图1-12library(Runuran)#为NIG产生分布目标distr1<-udghyp(lambda=-0.5,alpha=1,beta=0,delta=1,mu=0)gen1<-pinvd.new(distr1)x1<-ur(gen1,10000)dx1<-ud(gen1,x1)fx1<-up(gen1,x1)o1=order(x1)#为拉普拉斯分布产生分布目标distr2<-udlaplace(location=0,scale=1)gen2<-pinvd.new(distr2)x2<-ur(gen2,10000)dx2<-ud(gen2,x2)fx2<-up(gen2,x2)o2=order(x2)#为高斯分布产生分布目标x3<-rnorm(10000)dx3<-dnorm(x3)fx3<-pnorm(x3)o3=order(x3)#为柯西分布产生分布目标x4<-rcauchy(10000)dx4<-dcauchy(x4)fx4<-pcauchy(x4)o4=order(x4)#画图op<-par(mfrow=c(1,2))plot(sort(x1),dx1[o1],type="l",lty=1,xlab="X",ylab="Y",main="Distributioncomparison",xlim=c(-5,5),xaxp=c(-5,5,2),ylim=c(0,0.52),col="lightgreen")lines(sort(x2),dx2[o2],type="l",lty=2,col="black")lines(sort(x3),dx2[o3],type="l",lty=5,col="red")lines(sort(x4),dx2[o4],type="l",lty=3,col="blue")legend("topright",c("NIG","Laplace","Gaussian"),pch=c(1,1,1),col=c("lightgreen","black","red","blue"),cex=0.8,text.col=c("black","blue","red"),bty="n")plot(sort(x1),dx1[o1],type="l",lty=1,xlab="X",ylab="Y",main="Tailcomparation",xlim=c(-5,-4),xaxp=c(-5,-4,2),ylim=c(0,0.03),col="lightgreen")lines(sort(x2),dx2[o2],type="l",lty=2,col="black")abline(h=0,type="l",lty=5,col="red")lines(sort(x4),dx2[o4],type="l",lty=3,col="blue")points(x=-4.5,y=dcauchy(-4.5),cex=1,pch=1,col=1)text(x=-4.55,y=dcauchy(-4.55),"Cauchy",cex=0.8,col="blue",adj=c(0,-0.1))points(x=-4.5,y=ud(gen2,-4.5),cex=1,pch=1,col=1)text(x=-4.55,y=ud(gen2,-4.55),"Laplace",cex=0.8,col="black",adj=c(0,-0.1))points(x=-4.5,y=ud(gen1,-4.5),cex=1,pch=1,col=1)text(x=-4.55,y=ud(gen1,-4.55),"NIG",cex=0.8,col="lightgreen",adj=c(0,-0.1))points(x=-4.5,y=dnorm(-4.5),cex=1,pch=1,col=1)text(x=-4.55,y=dnorm(-4.55),"Gaussian",cex=0.8,col="red",adj=c(0,-0.1))par(op)1.8自助法例1.22x=rnorm(1000,0.1)y=pnorm(sort(x))plot(sort(x),y,type="l",col="red",xlab="X",ylab="edf(Y),cdf(X)",main="EDFandCDF,n=100")par(new=TRUE)f=ecdf(rnorm(100))plot(f,verticals=TRUE,do.points=FALSE,col="blue",main="",ylab="",axes=FALSE)图1-15col=c("blue","green")y=xx=x=matrix(1,100)x=rnorm(100)for(jin1:2){y=sample(1:100,100,replace=T)for(iin1:100){xx[i]=x[y[i]]}f=ecdf(xx)plot(f,verticals=TRUE,do.points=FALSE,col=col[j],main="",ylab="",axes=FALSE)par(new=TRUE)}f=ecdf(rnorm(100))plot(f,verticals=TRUE,do.points=FALSE,col="red",main="EDFand2bootstrapEDF's,n=100",ylab="edfs{1:3}(x)",xlab="x")

第2章多元正态分布理论2.4例2.3library("MASS")Sigma<-matrix(c(1,-1.2,-1.2,4),2,2)data=mvrnorm(n=200,c(3,2),Sigma)x=data[,1]y=data[,2]plot(x,y)

第3章基于因子的数据矩阵降维技术3.5实际应用#读入数据x1=read.table("G://data//carmark.txt",header=TRUE)x=as.matrix(x1)#数据标准化处理n=nrow(x)p=ncol(x)n1=rep(1,n)e=diag(1,n)h=e-1/n*n1%*%h%*%t(n1)s=1/n*t(x)%*%h%*%xns=nrow(s)d=diag(nrow=ns,ncol=ns)for(iin1:ns){d[i,j]=s[i,j]^(-1/2)}r=d%*%s%*%dxstar=1/sqrt(n)*h%*%x%*%d#求特征根xtx=t(xstar)%*%xstarlambda=eigen(xtx)$values#确定qpercentage=matrix(nrow=length(lambda),ncol=2)percentage[,1]=lambdapercentage[,1]=lambda/sum(lambda)#计算特征向量u1=eigen(xtx)$vectors[,1]u2=eigen(xtx)$vectors[,2]z1=xstar%*%u1z2=xstar%*%u2#画图plot(z1,z2,type="p",pch=18,col="maroon",cex=0.8,main="car")text(z1,z2,1:23,cex=0.8)abline(h=0,v=0)w1=sqrt(l1)*u1w2=sqrt(l2)*u2

第4章主成分分析#4.1library(graphics)z<-lm(dist~speed,data=cars)plot(cars)abline(a=80,b=-4)abline(z)abline(coef=coef(z))#例4.2library(stats)x=as.matrix(iris[1:100,1:4])xbar=colMeans(x)S=cov(x)ev=eigen(S)aa=eigen(S)$valuespca=princomp(x,cor=FALSE)y1=pca$scores[1:50,]y2=pca$scores[51:100,]par(mfrow=c(2,2))xmin=min(c(y1[,1],y2[,1]));ymin=min(c(y1[,2],y2[,2]))xmax=max(c(y1[,1],y2[,1]));ymax=max(c(y1[,2],y2[,2]))plot(y1[,1],y1[,2],xlab="pc1",ylab="pc2",xlim=c(xmin,xmax),ylim=c(ymin,ymax),main="第一个对第二个主成分",pch=19,col=3)par(new=TRUE)plot(y2[,1],y2[,2],xlab="pc1",ylab="pc2",xlim=c(xmin,xmax),ylim=c(ymin,ymax),main="第一个对第二个主成分",pch=17,col=2)xmin=min(c(y1[,1],y2[,1]));ymin=min(c(y1[,3],y2[,3]))xmax=max(c(y1[,1],y2[,1]));ymax=max(c(y1[,3],y2[,3]))plot(y1[,1],y1[,3],xlab="pc1",ylab="pc3",xlim=c(xmin,xmax),ylim=c(ymin,ymax),main="第一个对第三个主成分",pch=19,col=3)par(new=TRUE)plot(y2[,1],y2[,3],xlab="pc1",ylab="pc3",xlim=c(xmin,xmax),ylim=c(ymin,ymax),main="第一个对第二个主成分",pch=17,col=2)xmin=min(c(y1[,2],y2[,2]));ymin=min(c(y1[,3],y2[,3]))xmax=max(c(y1[,2],y2[,2]));ymax=max(c(y1[,3],y2[,3]))plot(y1[,2],y1[,3],xlab="pc2",ylab="pc3",xlim=c(xmin,xmax),ylim=c(ymin,ymax),main="第二个对第三个主成分",pch=19,col=3)par(new=TRUE)plot(y2[,2],y2[,3],xlab="pc2",ylab="pc3",xlim=c(xmin,xmax),ylim=c(ymin,ymax),main="第二个对第二个主成分",pch=17,col=2)plot(aa,xlab="指数",ylab="λ",main="S的特征根",pch=19)#4.3library(car)library(MASS)library(nnet)w=irisu=w[,-5]library(FactoMineR)res.pca=PCA(u,quanti.sup=3:4)#例4.5wine<-as.matrix(wine)[,-1]nwine<-scale(wine)R<-cor(nwine)GR<-svd(R)$ul<-svd(R)$dmodel<-princomp(nwine)summary(model)corr=GR%*%diag(l)^(1/2)#第一主成分及第二主成分的图解n=NULLplot(n,xlim=c(-1.1,1.1),ylim=c(-1.1,1.1),xlab="与第一主成分的相关系数",ylab="与第二主成分的相关系数")A=seq(0,2*pi,length.out=1000)r=1x1=r*cos(A)y1=r*sin(A)points(x1,y1,type="l")a=c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10","X11","X12","X13")for(iin1:13){text(corr[i,1],corr[i,2],a[i],cex=.7)}abline(h=0);abline(v=0)#178个样本的第一主成分和第二主成分散点图prin12<-nwine%*%GR[,1:2]plot(prin12[1:59,],xlim=c(-5,5),ylim=c(-4,4),col="blue",pch=2,xlab="第一主成分",ylab="第二主成分")points(prin12[60:130,],col="red",pch=3)points(prin12[130:178,],col="green",pch=16)

第5章因子分析#准备数据library("MASS")#Loadlibrary:Boston{MASS}bostonData=Boston[,-4]#Excludevariablex4bostonData=scale(bostonData)#Standardization#进行因子分析library("psych")#Loadlibrary:fa(psych),principal(psych)#MLM没有旋转mlmNonRotate=fa(bostonData,nfactors=3,fm="ml",rotate="none")#表5-3factor.plot(mlmNonRotate,labels=rownames(mlmNonRotate$loadings),show.points=FALSE)#图5-1#MLM具有varimax旋转mlmRotate=fa(bostonData,nfactors=3,fm="ml",rotate="varimax")#表5-4factor.plot(mlmRotate,labels=rownames(mlmRotate$loadings),show.points=FALSE)#图5-2#PCM具有varimax旋转pcmRotate=principal(bostonData,nfactors=3,rotate="varimax")#表5-5factor.plot(pcmRotate,labels=rownames(pcmRotate$loadings),show.points=FALSE)#图5-3#PfM具有varimax旋转pfmRotate=fa(bostonData,nfactors=3,fm="pa",rotate="varimax")#表5-6factor.plot(pfmRotate,labels=rownames(pfmRotate$loadings),show.points=FALSE)#图5-4

第6章聚类分析#例6.1x1=c(0,0)x2=c(1,0)x3=c(6,8)X=rbind(x1,x2,x3)dist(X,method="manhattan",upper=TRUE)A=dist(X,upper=TRUE)#计算欧式距离D=A*A#平方的欧式范数#例6.2X=USArrests[1:12,]A=dist(X,upper=TRUE)D=A*A#6.4data(iris)attach(iris)iris.hc<-hclust(dist(iris[,1:4]),"complete")plclust(iris.hc,labels=FALSE,hang=-1)iris.id<-cutree(iris.hc,3)table(iris.id,Species)s1=c(0,0,0,0)s2=c(0,0,0,0)s3=c(0,0,0,0)for(iin1:150)s1=iris[i,1:4]*(iris.id[i]==1)+s1for(iin1:150)s2=iris[i,1:4]*(iris.id[i]==2)+s2for(iin1:150)s3=iris[i,1:4]*(iris.id[i]==3)+s3s1=s1/50s2=s2/72s3=s3/28iris.hc<-hclust(dist(iris[,1:4]),"average")plclust(iris.hc,labels=FALSE,hang=-1)iris.hc<-hclust(dist(iris[,1:4]),"ward")plclust(iris.hc,labels=FALSE,hang=-1)

第7章判别分析7.1.1极大似然判别准则例7.3x<-seq(-3,3,0.01)y1<-dnorm(x,0,1)y2<-dnorm(x,1,1/2)v1=1/3*(4-sqrt(4+6*log(2)))v2=1/3*(4+sqrt(4+6*log(2)))plot(x,y1,type="l",ylim=c(0,0.8),ylab="密度函数")lines(x,y2)abline(v=v1,lty=2)abline(v=v2,lty=2)text(-1.5,0.6,"R1",col=2)text(2.8,0.6,"R1",col=2)text(1.7,0.7,"R2",col=4)7.2.1错判概率的估计例7.4x<-as.matrix(iris[1:4])x1<-x[1:50,];x2<-x[51:100,];x3<-x[101:150,]xbar1<-apply(x1,2,mean)xbar2<-apply(x2,2,mean)xbar3<-apply(x3,2,mean)su<-50/(150-3)*(var(x1)+var(x2)+var(x3))h12<-function(x)t(xbar1-xbar2)%*%solve(su)%*%(x-1/2*(xbar1+xbar2))h13<-function(x)t(xbar1-xbar3)%*%solve(su)%*%(x-1/2*(xbar1+xbar3))h23<-function(x)t(xbar2-xbar3)%*%solve(su)%*%(x-1/2*(xbar2+xbar3))actual<-c(rep(1,50),rep(2,50),rep(3,50))predicted<-NULLfor(iin1:150){if(h12(x[i,])>=0&&h13(x[i,])>=0)predicted[i]=1elseif(h12(x[i,])<0&&h23(x[i,])>=0)predicted[i]=2elsepredicted[i]=3}table(actual,predicted)

第8章对应分析#例8.2#多元对应分析-例8.2比利时人读报数据newspaper<-read.table("C:/Users/TIAN/Desktop/newspaper.txt",header=TRUE)#多元对应分析-例8.2比利时人读报数据newspaper