混料试验的正交格点填充设计①

李光辉， 李俊鹏，， 张崇岐

1.凯里学院 理学院，贵州 凯里 556011； 2.广州大学 经济与统计学院，广州 510006

混料试验设计是一类特殊的设计， 它广泛应用于生物、 医学、工业以及科研实践中[1-4]．混料试验的试验域是q-1维正规单纯形

为了使得试验点能均匀分布在试验域内， 许多统计学者尝试了各种不同的方法．文献 [5]介绍了混料试验设计的基本概念及应用．关于混料均匀设计方面的研究可参见文献[6-9]．

文献[10-11]研究了混料试验域上构造均匀设计的方法； 文献 [12]讨论了在具有附加约束的混料试验域内构造均匀设计的方法； 文献[13]给出了一种单纯形内的切片设计． 这些方法所研究的设计都是在单纯形内部分布均匀， 且在某种均匀性准则下能达到较好的水准． 由于混料试验通常使用回归模型来建模， 所以均匀设计具有很好的稳健性．

有时需要在混料试验域内采用非参数回归模型来建模， 除了需要大量的试验点外， 还需要满足“试验点之间的距离尽可能相等”这一条件． 混料格子点集虽然能均匀分布在混料试验域的内部及边界， 但是当混料分量数超过3以后， 相邻的格点之间的距离是不相等的， 并且由混料格子点集剖分得到的子单纯形也不完全相等． 这使得使用混料格子点集构造的设计不利于使用非参数回归模型来建模． 参考超立方体内的建模方法， 正交格子点集不仅排布均匀， 相邻的格点之间的距离相等， 并且超立方体内的完全析因设计有诸多优良的性质．

关于超立方体上的格点设计已有一些研究成果． 文献[14-15]研究了交错网格的极大极小距离设计， 以及旋转球填充设计．这些空间填充设计都是基于格子点集来研究的， 并且在超立方体上有很多好的性质．

鉴于此， 我们考虑将超立方体内的正交格子点集变换到混料试验域内， 并且保证变换后试验点之间的距离不变． 根据正交格子点集的性质， 研究在混料试验域上正交格点设计的性质．首先给出一种线性变换， 能将Rq-1的点集变换到Sq-1中， 并且保证点之间的距离不变； 然后介绍在单纯形内部， 将试验点以旋转或添加的方式， 使得设计在混料试验域内更为均匀； 最后给出实例分析， 说明单纯形内的正交格点设计具有很好的均匀性， 并且非参数建模也具有稳健性．

1 基本记号

对于q分量混料系统， 若试验还受到其他约束条件的限制， 我们将试验域记作

X={x=(x1，x2，…，xq)T：x∈Sq-1，aj≤φj(x)≤bj，j=1，2，…，k}

其中：φj(x)是关于x∈Sq-1的已知函数；aj，bj(j=1，2，…，k)为已知常数．通常X是Sq-1内部的一个不规则凸几何体．

对于任意的x=(x1，x2，…，xq)T∈Sq-1⊂X以及一个n点设计

Pn={x1，x2，…，xn}⊂X

我们给出点到点集的距离定义为

其中d2(x，xi)=‖x-xi‖2为欧氏距离的平方．

文献[16]提出了3种度量混料试验域设计的均匀性指标， 即偏差的定义． 为方便讨论， 我们主要讨论均方误偏差和最大距离偏差．

1)均方误偏差

2)最大距离偏差

以上两种偏差的计算较为复杂， 在实际情形下一般都是采用Monte Carlo方法来计算近似值． 首先生成试验域内的随机混料点集， 或使用高阶格子点集

令LX=L{q，m}∩X={t1，t2，…，tN}作为度量均匀性的NT-net， 再令

(1)

使用msed(Pn)和md(Pn)分别作为MSED(Pn)和MD(Pn)的近似值．在实际计算中， 一般要求N>1 000．

q-1维空间中的正交格子点集可以表示为

(2)

其中：αj1，αj2，…，αj，q-1∈Z；eq-1(i)为q-1维向量， 其第i个元素为1， 其余元素为0；c是已知的正数， 它表示正交格点之间的距离． 我们记由(2)式生成的正交格子点集为

对于任意的y∈H(q-1，c)， 定义y0的voronoi单元为

Vor(y0)={z：‖z-y0‖2≤‖z-y‖2，∀y∈H(q-1，c)}

这样， 我们可以考虑将超立方体内的试验点集转化到单纯形中， 并且保持格点之间的距离不变．

2 保距独立变换

为了实现将q-1维超立方体内的正交格子点集变换到单纯形中， 我们构造一种线性变换， 将Rq-1中的点变换到Sq-1中， 且点与点之间的距离保持不变．文献[17]的单纯形保距独立变换方法如下：

定义一个q-1阶矩阵

(3)

WTW=Iq-1+Jq-1

对于W的任意两列wi，wj，i≠j， 有

这种变换可以将单纯形Sq-1中的点映射到Rq-1中．

对于任意的xi=(xi1，xi2，…，xiq)T∈Sq-1， 记x-1，i=(xi2，xi3，…，xiq)T， 即yi=Wx-1，i，i=1，2， 有‖y1-y2‖=‖x1-x2‖恒成立．

这是因为由矩阵W的性质可知

由此可见， 经过变换后任意两点在Rq-1中的距离与Sq-1中是一致的．

3 填充设计

我们按照以下方式生成Sq-1中的正交格点．

首先生成Rq-1中的正交格子点集{y1，y2，…，yN}⊂H{q-1，c}， 在混料试验域中选择一点x0=(x01，x02，…，x0q)T∈X作为参照点， 令y0=Wx-1，0，其中

x-1，0=(x02，x03，…，x0q)T

其次， 令矩阵

Y=(y1，y2，…，yN)T+1Ny0

(4)

是以y0作为原点而形成的正交格子点集．

最后，进行线性变换， 令矩阵

X=(1N-X-11q-1，X-1)=(x1，x2，…，xN)T

(5)

HX=H ′{q-1，c}∩X

(6)

表示最终确定落入试验域内的正交格子点集．

(x1，x2，…，x2q)T=((a1，b1)T⊗12q-1，12⊗(a2，b2)T⊗12q-3，…，12q-1⊗(aq，bq)T)

则有以下定理．

证首先证明对于任意点x0=(x01，x02，…，x0q)T∈H0， 有

为了使得试验域X内的正交格子点集具有良好的均匀性， 即使得设计

Pn={x1，x2，…，xn}=H ′{q-1，c}∩X

的均方误偏差MSED(Pn)与最大距离偏差MD(Pn)尽可能地小， 我们考虑从以下两个方面来实现．

4 旋转变换

在X中实施旋转变换， 使得落入X内的正交格点尽可能地多． 我们令q-1阶方阵为

Rij(θij)=(eq-1(1)，…，eq-1(i-1)，Ei，eq-1(i+1)，…，eq-1(j-1)，Ej，eq-1(j+1)，…，eq-1(q-1))

其中Ei=eq-1(i)cosθij-eq-1(j)sinθij，Ej=eq-1(j)sinθij+eq-1(j)cosθij．在(4)式的基础上， 再令

其中Θ=(θ12，θ13，…，θ(q-2)(q-1))为旋转参数向量． 将Y(Θ)经过(5)式的逆变换后得

(x1，x2，…，xN)=W-1YT(Θ，c)

(7)

令HX={x1，x2，…，xN}∩X表示最终确定落入试验域内的正交格子点集． 通过调整(7)式中的参数Θ， 使得HX达到满足要求的均匀性．

5 补充设计点

首先， 在X内使用高阶混料格子点集

LX=L{q，m}∩X={t1，t2，…，tN1}

作为在X上度量一个设计所使用的NT-net． 令

这里uij(i=1，2，…，N；j=1，2，…，q)相互独立， 且都服从均匀分布uij～U(0，1)．

使用文献[19]的拟分量变换方法， 选择适当的拟分量变换参数λ>0， 令

msed{Pn，xn+1}→inf

Pn+k={x1，x2，…，xn，…，xn+k}

(8)

考虑在试验域X上使用正交格子点集来建模， 按(7)式或(8)式生成X上的正交格子点集记为PN={x1，x2，…，xN}．定义高斯核函数为

其中h为已知的窗宽参数， 若试验在各点处的观测值为y1，y2，…，yN， 则基于正交格点的非参数回归模型为

如果在试验域内有足够多的格点， 选择适合的窗宽， 能使得非参数回归模型具有很好的预测效果．

例1在S3-1上构造一个正交格点设计， 使得各相邻点之间的距离为0.2， 通过旋转过后设计点的分布以及试验域内各个点到该设计的距离等高线图如图1所示．

图1 S3-1内旋转后的正交格点设计

以模型

为例， 在S3-1上的曲面图如图2(a)所示， 使用正交格点建立非参数回归模型， 拟合回归曲面如图2(b)所示．

图2 曲面图

由图2可见， 使用正交格点设计建立的非参数回归模型具有很好的稳健性， 拟合效果还是不错的．

使用正交格点构造的混料试验域内的设计分布均匀，利于非参数建模， 使用文献[20]的方法， 结合正交格子点集能有效检验设计的最优性， 关于这方面的性质有待进一步研究．