多元统计，实例分析(置信语言样本区间函数)

library(mvtnorm)

Sigma<-matrix(c(10,3,3,2),ncol=2)

set.seed(10)

x<-rmvnorm(n=30000,mean=c(1,2),sigma=Sigma) #生成样

本容量 30000 的二维正态样本

plot(x,pch=20,cex=0.7,xlab="X1",ylab="X2") #基本散点图，

难以反映点的密度

smoothScatter(x,xlab="X1",ylab="X2") #生成密度估计图

#R语言#

R语言还提供了三维图形的作图函数，如“rgl”数据包的“plot3d”函数，该函数能创建交互式的三维散点图，使用者可通过鼠标操作，旋转、放缩三维图形。
可以利用“dmvnorm”函数计算样本各点的概率密度值，以“plot3d”函数在三维图中标示出每一点的概率密度。
采用图1的数据，制作其概率密度的交互式三维散点图。

library(rgl)

z=dmvnorm(x,mean=c(1,2),sigma=Sigma) #计算概率密度

plot3d(x[,1],x[,2],z,xlab="X1",ylab="X2",zlab="density",cex=0.8)

根据“plot3d”函数生成的交互式三维图的两个旋转方向，先将三维图旋转至X1-X2平面，此时观察到的样本点二维分布与图1效果一致，再将三维图旋转，看见二维散点图转换成三维概率密度图的动态过程，形成比较深刻的视觉效果。

#计算机辅助化工装置选型设计#

利用R语言“mvtnorm”包中的“rmvnorm（）”函数生成二维正态样本，再利用“plotrix”数据包中的“draw.ellipse（）”函数制作二维置信域（置信椭圆）。
“draw.ellipse（）”函数的主要参数有：draw.ellipse(x,y,a=1,b=1,angle=0,deg=TRUE,...)其中，x、y为椭圆中心，还有样本均值向量的坐标；a为长轴半径，λ1是样本协差阵S的最大特征值；b为短轴半径，λ2是样本协差阵S的最小特征值；angle是椭圆长轴与横轴夹角，即λ1对应的单位特征向量与横轴的夹角，利用反正切函数即可得到该夹角取值；deg取值为TRUE时，表示该夹角为角度制，取值FALSE，夹角为弧度制。
相关代码如下所示。

library(plotrix)#载入plotrix包

library(mvtnorm)#载入mvtnorm包

Sigma<-matrix(c(3,2,2,2),2)#生成总体协方差阵

set.seed(5)

x<-rmvnorm(n=30,c(0,1),Sigma)#生成样本容量为30的二维正态随机样本

xbar<-colMeans(x)#计算样本均值向量

S<-cov(x)#计算样本协方差阵S

eigenX<-eigen(S)#计算S的特征值与单位特征向量

lamda1<-eigenX$values[1]#读取最大特征值

lamda2<-eigenX$values[2]#读取第二大特征值

n<-nrow(x)#读取样本容量

p<-ncol(x)#读取维数

Fa<-qf(0.95,p,n-p)#计算相应的F分布的上分位数

a1<-sqrt(Fap(n-1)/(n-p)lamda1/n)#计算长轴半径

b1<-sqrt(Fap(n-1)/(n-p)lamda2/n)#计算短轴半径

angle=atan(eigenX$vectors[2,1]/eigenX$vectors[1,1])#通过反正切函数计算夹角

plot(x,pch=19,cex=0.6,xlim=c(-3.5,4.5),ylim=c(-3.5,4.5),xlab="X1",ylab="X2")#画二维散点图

draw.ellipse(xbar[1],xbar[2],a1,b1,angle,deg=FALSE,lwd=1.8)#编制置信椭圆

legend("bottomright","95%置信域",lty=7,cex=0.7)#编制图例

图5中，椭圆区域即为该组数据的总体期望μ的95%置信域。
此外，根据式（3）和式（4），可以在图5中增添相应的95%置信度的T2联合置信区间和邦弗伦尼联合置信区间，由于联合置信区间是一个矩形，因此可利用“rect（）”函数进行添加，相关代码如下所示。

#T2联合置信区间

Fa<- qf(0.95,p,n-p)#计算相应的 F 分布的上分位数

td1<-sqrt(Fap(n-1)/(n-p))sqrt(S[1,1]/n) #第一个分量的

置信区间宽度

td2<-sqrt(Fap(n-1)/(n-p))sqrt(S[2,2]/n) #第二个分量的

置信区间宽度

tx1<-xbar[1]-td1 #第一个分量的置信下限

tx2<-xbar[1]+td1 #第一个分量的置信上限

ty1<-xbar[2]-td2 #第二个分量的置信下限

ty2<-xbar[2]+td2 #第二个分量的置信上限

rect(xleft=tx1, ybottom =ty1, xright =tx2, ytop =ty2, lty=5)

#制作 T

2联合置信区间

#邦弗伦尼联合置信区间

bd1<-qt(1-0.05/4,n-1)sqrt(S[1,1]/n) #第一个分量的置信

区间宽度

bd2<-qt(1-0.05/4,n-1)sqrt(S[2,2]/n) #第二个分量的置信

区间宽度

bx1<-xbar[1]-bd1 #第一个分量的置信下限

bx2<-xbar[1]+bd1 #第一个分量的置信上限

by1<-xbar[2]-bd2 #第二个分量的置信下限

by2<-xbar[2]+bd2 #第二个分量的置信上限

rect(xleft=bx1, ybottom =by1, xright =bx2, ytop =by2, lty=3)

#制作邦弗伦尼联合置信区间

legend("bottomright",c("95%置信域","95%联合 T2 置信区

间","95%邦弗伦尼联合置信区间"),lty=c(7,5,3),cex=0.7) #编制图例

代码的编程结果如图6所示。
直观了解二元正态分布期望μ的联合T2置信区间是其联合置信域的外接矩形，由于联合置信域的置信度是95%，因此联合T2置信区间的置信度将明显高于95%;直观了解邦弗伦尼联合置信区间在保证置信度的基础上，区域范围一般小于联合T2置信区间；联合置信域中椭圆长短轴的半径、方向的计算方法；熟悉了R语言中“mvrnorm（）”“eigen（）”“draw.ellipse（）”“rect（）”等函数的使用方法。