很早就想为COS写一篇关于计算机试验的东西。可是始终也未敢动笔,觉得自己才疏学浅,生怕写得偏颇。但是另外一方面,又觉得这是COS上一块空白的话题,没有人提及过。今天写这篇文章,主要是为抛砖引玉,另外丰富一下主站文章的题材。

在这篇文章中,我想大致介绍一下计算机试验的设计以及建模,另外会有一些R中专门做计算机试验的包的相关介绍。

目前,比较流行的计算机试验设计与建模的教材有两本(似乎也只有这两本)。一本是Thomas J. Santner, Brian J.Williams以及William I.Notz合著的,Springer2003年出版的《The Design and Analysis of Computer Experiment》。另外一本是Kai-Tai Fang ,Runze Li以及Agus Sudjianto合著的,2005年由Chapman&Hall出版的《Design and Modeling for Computer Experiment》。所谓的计算机试验(computer experiment),是相对于传统的实体试验(physical experiment)而言的。我们都知道,实体试验的数据(比如农业试验,生物试验等得到的数据)总是会受到随机误差的影响,因此在实体试验设计中,引入了“重复”的原则,目的就是要通过重复的试验,减小随机误差对于分析的影响。(关于实体试验的设计以及建模,国内外已经有很多教材讨论这些,我相信大多数统计专业的同学所学习的试验设计入门课程都是在学习实体试验的设计以及数据分析。)与此相对的,计算机试验得到的数据并不受到随机误差的干扰。因为一段固定的代码,在计算机上无论运行多少次,得到的结果都是一样的。

计算机试验,通常比实体试验包含了更多的变量,用于研究特别复杂的系统,比如航天探测器。对于这些特别复杂的系统,往往需要用一个更加简单的拟模型(meta model)来逼近。计算机试验的目的之一,就是要通过在计算机上,利用代码模拟某个复杂的设备或者某个过程,得到了试验的数据之后,通过分析建模,从而得到一个相对更简单的模型,被称为拟模型,拟模型是非常有用的。

由于计算机试验具有复杂性,以及输出结果的确定性,做计算机试验收集数据的时候,我们需要一些特殊的设计。

试验设计,总是和模型紧密相联系的。事实上,试验设计是一种收集数据的策略。在实体试验中,当我们需要比较某种因子不同水平之间的差别,我们通常会选用因析设计。因析设计事实上就是和因子模型相对应而产生的收集数据的策略。另外,所谓的最优设计,也是相对应于模型的最优化准则而产生的收集数据的策略。而在计算机试验中,常用的是模型是所谓的全局均值模型(overall mean model)。

不失一般性,我们考虑试验区域是一个s维的单位立方\(C^s=[0,1]^s\)。当试验次数给定了,比如说\(n\)次,我们就需要考虑如何找到一个好的试验\(D\_n=\{x\_1,\cdots,x\_n\}\),使得\(f(x)-g(x)\)对于试验区域上的所有点,都尽可能的小。其中\(f(x)\)是真实的模型,而\(g(x)\)是拟模型。我们考虑到,对于所有的点,都要求上述差分达到最小是很困难的。因此,很多学者都将该问题进行了简化,转而去寻找全局均值\(E(y)=\int\_{C^s}f(x)dx\)的最佳的估计值。而通常,都用样本均值,也就是平均数\(\bar{y}(D\_n)=\frac{1}{n}\sum\_{i=1}^nf(x_i)\)来估计总体均值。因此,我们的问题也就转化为了如何找到一个试验,使得这个试验得到的数据做样本均值来估计总体均值是最好的。这也就是计算机试验设计的动因。

从统计学的角度来看,如果样本点\(x\_1,\cdots,x\_n\)是在\(C^s\)上均匀分布上独立抽取的,那么样本均值是无偏的,并且方差为\(\frac{var(f(x))}{n}\)。但是,通常来说,这个方差太大了。因此很多学者提出了不同的抽样方法来降低这个样本均值的方差。1979年,McKay,Beckman和Conover基于分层抽样提出了一种新的抽样方法,这种方法使得随机选择的样本点并不独立,并且边际分布相同,因此

$$ Var(\bar{y}(D_n))=\frac{1}{n}Var(f(x))+\frac{n-1}{n}Cov(f(x_1),f(x_2)) $$

\(f(x\_1)\)\(f(x\_2)\)负相关时,方差就减小了。这种方法被人不断的推广,创新。成为了现在计算机试验设计中一块非常重要的内容——拉丁超立方体设计(Latin-Hypercube design)。在后来的不断发展中,拉丁超立方和正交表相结合,产生了很多正交拉丁超立方设计。此外,当引入一些准则时,比如熵准则,最大最小距离准则等,也产生了对应该准则的最优的设计,包括最大最小拉丁超立方等。这些设计在计算机试验中应用非常广泛,最近也有人将这种设计用在变量选择的过程中,通过模拟,大大提高了变量选择的精确度。

与拉丁超立方设计相对应的另外一种设计,是由我国科学家方开泰,王元在上世纪八十年代提出的均匀设计。均匀设计这一想法的提出,起源于伪蒙特卡洛方法中的Koksma-Hlawaka不等式:\(|E(y)-\bar{y}(D\_n)|\leq V(f)D(D\_n)\)。其中\(D(D\_n)\)是设计\(D\_n\)的星偏差(star-discrepancy),而\(V(f)\)是方程\(f\)的总体变动情况。其中星偏差是一种均匀性的度量,当试验点在试验区域分布越均匀时,星偏差的值就越小,因此,样本均值与总体均值之差的绝对值的上界也就越小,相当于提高了估计的精度。一个试验区域上具有最小星偏差的设计,就是所谓的均匀设计。均匀设计由我国科学家自主提出,并且在国际上得到了非常好的反响以及应用。之后,有很多学者,对于均匀设计做出了改进,使得均匀设计也不断地发展壮大。

下面简单介绍一下对于计算机试验得到的数据建模的常用方法。事实上,统计学习中的常用的建模方法都可以用来对计算机试验数据进行建模。包括样条回归,局部回归,神经网络等,都可以对数据进行建模。这里要要重点介绍的是一种非常常用的模型,高斯-克里金模型(Gaussian-Kriging models)。克里金是南非地质学家,他在其硕士论文中分析矿产数据的时候,提出了克里金模型,后来由几位统计学家发展了他的模型,得到了目前在计算机试验中非常常用的高斯-克里金模型。

$$ y(x)=\sum_{j=0}^{L}\beta_jB_j(x)+z(x) $$

其中\(B_j(x)\)是选定的基函数,\(z(x)\)是随机误差,来自于高斯过程。该高斯过程零均值,方差为\(\sigma^2\),协方差函数为\(r(\theta,s,t)=corr(z(s),z(t))\)。通常情况下\(r\)是预先选定的。可以证明,高斯克里金模型的预测值是BLUP。

关于计算机试验数据的建模,方法很多,但是建模方法仅仅是建模方法,无论是计算机试验还是实体试验,建模方法往往通用,因此,这里不做展开论述了。

下面简单介绍一下R语言中,做计算机试验的一些包。

lhs包,主要用来生成拉丁超立方体设计。DiceDesign包可以生成更多的一些空间填充设计。DiceKriging提供了对计算机试验数据构建高斯-克里金模型的功能,DiceView则提供了多维模型可视化的手段。以上的包在task view中都有所介绍。此外,值得一提的是VizCompX包,其对计算机试验的数据拟合高斯模型,并且提供可视化的手段。感兴趣的童鞋可以参看这几个包的帮助文档,帮助文档都不长,容易上手。

最后,提一点我个人的小感想。我认为,计算机试验中,更应该引人注意的话题应该是试验的设计。一个好的设计,不光能够提高试验的精度,更能体现一种美感。生成一个好的设计,大致上有两种路可以走,一种是预先设定一个准则,然后利用最优化的算法,进行搜索。常用的算法包括一些全局的最优算法,比如门限接受算法,模拟退火算法等。另外一种方法,则是利用数学中其他学科,比如伽罗瓦理论,编码理论等,通过这些学科的理论特点,精巧地构造一个设计。这两种构造设计的方法相互补充,相辅相成,提供了非常广大的研究空间。另外,这些为计算机试验而提出的设计,也可以用在很多其他的统计学的课题中,比如之前已经提到过的变量选择方法精确度的提升方面。总之,试验设计是研究如何有效收集数据的,无论数据是来自于实体试验,还是计算机试验,一个设计了的试验,往往都可以提高效率,提高精度。这也是试验设计这门课程的意义所在。越南的Nyuan开发了一款试验设计的软件Gendex,通过输入试验的行数,列数,选择你所需要的设计种类,便可以生成多种不同的设计阵。这款软件,是目前可以生成设计阵种类最多的,感兴趣的童鞋,不妨搜来玩一玩。

发表/查看评论