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

1、多元统计分析代码第1章 多元分布1.5 样本分布和极限定理例1.8 n=5 #n=50X=matrix(0,n,1000)for (i in 1:1000)X,i=rbinom(n,1,1/2)XX=(colMeans(X)-1/2)*5(1/2)plot(density(XX),type=l,xlab=1000 Random.Sample, ylab=Estimated and normal Density,main=Asymptotic Distribution, N=5,col=blue,xlim=c(-2,2)par(new=TRUE)dat=rnorm(100000,0,1/2)pl

2、ot(density(dat),col=red,main=,axes=FALSE)1.6 厚尾分布curve(dnorm(x),-5,5,xlab=X,ylab=Y,main=distribution comparison, 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(b

3、lue,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)#制造生成器目标；利

4、用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),dx

5、1o1,type=l,xlab=X,ylab=Y, main=pdf of GH,HYP and NIC,xlim=c(-5,5), xaxp=c(-5,5,2),ylim=c(0,0.52),col=black)lines(sort(x2),dx2o2,type=l,col=red)lines(sort(x3),dx2o3,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

6、o1,type=l,xlab=X,ylab=Y, main=cdf of GH,HYP and NIC,xlim=c(-5,5),xaxp= c(-5,5,2),ylim=c(0,1),col=black)lines(sort(x2),fx2o2,type=l,col=red)lines(sort(x3),fx2o3,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-5

7、op-par(mfrow=c(1,2)curve(dt(x,3),-5,5,xlab=X,ylab=Y,main=pdf of t-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),c

8、ol=c(black,blue,red) ,text.col=c(black,blue,red),bty=n)curve(pt(x,3),-5,5,xlab=X,ylab=Y,main=cdf of t-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

9、,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=tail comparison-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)curv

10、e(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=deep pink)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,c

11、ol=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=deep pink)text(x=3.21,y=dt(3.2,4

12、5),t45,cex=0.8,col=deep pink,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=tail comparison-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=

13、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

14、(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)#产生大小为100

15、00的样本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(ge

16、n3,x3)o3=order(x3)#绘图op-par(mfrow=c(1,2)plot(sort(x1),dx1o1,type=l,xlab=X,ylab=Y, main=pdf of Laplace distriction,xlim=c(-5,5), xaxp=c(-5,5,2),ylim=c(0,0.52),col=black)lines(sort(x2),dx2o2,type=l,col=blue)lines(sort(x3),dx2o3,type=l,col=red)legend(topright,c(L1,L1.5,L2),pch=c(16,16,16), col=c(black,

17、blue,red),text.col=c(black,blue, red),bty=n)plot(sort(x1),fx1o1,type=l,xlab=X,ylab=Y, main=cdf of Laplace distribution,xlim=c(-5,6),xaxp= c(-5,5,2),ylim=c(0,1),col=black,pch=20)lines(sort(x2),fx2o2,type=l,col=blue)lines(sort(x3),fx2o3,type=l,col=red)legend(topleft,c(L1,L1.5,L2),pch=c(16,16,16),col=

18、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=pdf of Cauchy distribution, 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

19、,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=cdf of Cauchy distribution, xaxp=c(-5,5,2),ylim=c(0.09,0.91),col=black)curve(pcauchy(x,0,1.5),-5,5,add=

20、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,

21、xlabb=X,ylab=Y, main=pdf of Gaussian mixture and Gaussian,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(Gaussian mixture,Gaussian),pch=c(16,16), col=c(blue,red), text.col=c(blue,red)

22、,bty=n)curve(0.8*pnorm(x)+0.2*pnorm(x,0,3),-9,9,xlabb=X,ylab=Y, main=cdf of Gaussian mixture and Gaussian,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(Gaussian mixture,Gaussian),pch=c(16,16),

23、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产

24、生分布目标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产生

25、分布目标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),dx1o1,type=l,xlab=X,ylab=Y, main=tail comparison-GH,xlim=c(4,5), xaxp=c(4,5,2),ylim=c(0,0.02),col=black)lines(sort(x2),dx2o2,

26、type=l,col=red)lines(sort(x3),dx2o3,type=l,col=blue)lines(sort(x4),dx2o4,type=l,col=deep pink)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),HY

27、P,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=deep pink)text(x=4.55,y=ud(gen4,4.55),NIG,cex=0.8, col=deep pink,adj=c(0,-0.1)plot(sort(x1),(x10.5*exp

28、(-x1)o1,type=l,xlab=X,ylab=Y, main=tail comparision-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=deep pink)points(x=4.5,y=4.50.5*exp(-4.5)

29、,cex=0.8,pch=1,col=black)text(x=4.55,y=4.550.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)

30、,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=deep pink)text(x=4.55,y=4.55(-1.5)*exp(-4.55),NIG,cex=0.8,col=deep pink,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-u

31、r(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)f

32、x4-pcauchy(x4)o4=order(x4)#画图op-par(mfrow=c(1,2)plot(sort(x1),dx1o1,type=l,lty=1,xlab=X,ylab=Y, main=Distribution comparison,xlim=c(-5,5), xaxp=c(-5,5,2),ylim=c(0,0.52),col=light green)lines(sort(x2),dx2o2,type=l,lty=2,col=black)lines(sort(x3),dx2o3,type=l,lty=5,col=red)lines(sort(x4),dx2o4,type=l,l

33、ty=3,col=blue)legend(topright,c(NIG,Laplace,Gaussian),pch=c(1,1,1),col=c(light green,black,red,blue) ,cex=0.8,text.col=c(black,blue,red),bty=n)plot(sort(x1),dx1o1,type=l,lty=1,xlab=X,ylab=Y, main=Tail comparation,xlim=c(-5,-4), xaxp=c(-5,-4,2),ylim=c(0,0.03),col=light green)lines(sort(x2),dx2o2,type

34、=l,lty=2,col=black)abline(h=0,type=l,lty=5,col=red)lines(sort(x4),dx2o4,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=

35、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=light green,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.22 x=rnorm(1000,0.1)y=

36、pnorm(sort(x)plot(sort(x),y,type=l,col=red,xlab=X,ylab=edf(Y),cdf(X), main=EDF and CDF,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(j in 1:2)y=sample(1:100,100,replace=T)for(i in 1

37、:100)xxi=xyif=ecdf(xx)plot(f,verticals =TRUE,do.points=FALSE,col=colj, main=,ylab=,axes=FALSE)par(new=TRUE)f=ecdf(rnorm(100)plot(f,verticals = TRUE,do.points=FALSE,col=red, main=EDF and 2 bootstrap EDFs,n=100,ylab=edfs1:3(x),xlab=x)第2章 多元正态分布理论2.4例2.3library(MASS)Sigma - matrix(c(1,-1.2,-1.2,4),2,2)

38、data=mvrnorm(n=200, c(3, 2), Sigma)x=data,1y=data,2plot(x,y)第3章 基于因子的数据矩阵降维技术3.5 实际应用#读入数据x1=read.table(G:/data/car mark.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(i in 1:ns) di,j=si,j(-1

39、/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,1u2=eigen(xtx)$vectors,2z1=xstar%*%u1z2=xstar%*%u2#画图plot(z1,z2,type=p,pch=18,col=ma

40、roon,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(iris1:100,1:4)xbar=colMeans(x)S=cov(x)ev=eigen(S)aa=eigen(S)$v

41、aluespca=princomp(x,cor=FALSE)y1=pca$scores1:50,y2=pca$scores51: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,xla

42、b=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);y