在考虑多个变量存在内部联系的条件下，挖掘、分析多维总体内在的数量特征与统计规律。R语言是属于GNU系统的一个自由、免费、源代码开放的软件，具有一套完整的数据处理、计算和制图功能。近年来，R语言已逐步成为统计工作者一种重要的分析软件。R语言在多元统计分析的几个优势：（1）R语言是免费软件。R语言是一种完全免费且源代码开放的软件，使用者可以在R语言的官方网站或者其它资源网站中下载安装程序、数据包和相关的文档材料。相比于其他的商业统计软件，更容易在自己电脑上安装R语言。（2）R语言是可编程的软件。作为一个开放的统计编程软件，R语言允许使用者自行编程，而且R的语法简单灵活，容易掌握。熟悉R的语法后，便可根据需要，自行编制R语言的例子。(3)R语言更新速度快。基于R语言的可编程性，使用者可以编制自己的函数来扩展现有的分析算法，因此R语言可实现的统计算法更新速度明显高于SPSS、SAS等软件，大多数最新的统计方法和技术都可以在R中直接实现。（4）R语言的可视化功能强大。R语言提供了大量的绘图函数，使用者只需调用合适的绘图函数，并对绘图参数进行相应设置，便可得到各式各样的图形结果，无论是基本的直方图、饼图、条形图等，还是复杂的组合图、三维图、动画等，都可以采用R语言实现。可以利用R语言的“mvtnorm”数据包中的“rmvnorm（）”函数，生成一组样本容量为30000的二维正态分布随机样本，并作出该样本的散点图。服从二维正态分布的随机样本散点图呈现出椭圆的形状，但是否“越靠近中心，点的密度越大”，这一点仍需继续验证。该问题可利用“smoothScatter（）”函数生成该组样本的密度估计图，越靠近样本中心的颜色越深，说明该区域的点密度越高，而且相同密度的区域边界呈现出明显的椭圆形状，由此可以把概率密度直观展示出来。相关代码如下所示。

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(Fa*p*(n-1)/(n-p)*lamda1/n)#计算长轴半径

b1<-sqrt(Fa*p*(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(Fa*p*(n-1)/(n-p))*sqrt(S[1,1]/n) #第一个分量的

置信区间宽度

td2<-sqrt(Fa*p*(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（）”等函数的使用方法。