方开泰、刘民千、周永道《试验设计与建模》课件-5

均匀设计1目录引言总体均值模型均匀性度量-偏差均匀设计的构造好格子点法及其推广随机优化法均匀设计的应用正交性与均匀性的联系25.1引言传统试验设计中的未知参数正交设计：主效应和二阶交互效应的个数为最优回归设计：多元二次模型的未知参数个数为3模型未知设计：均匀布点建模：寻找近似模型g:未知函数;

x1,

,xs

:s

个因素;T=[a1,b1]

[as,bs]:试验区域,

其中

aj

xj

bj,j=1,

,s4例1.6(续)事先假设模型未知试验次数n=12,取q(=2,3,4,6,12)个不同的均匀的试验点，并作适当重复(n/q

次)。采用d阶多项式回归模型重复N=501次试验，得到各模型的N个MSE。取中位数对应的模型如下图565.2总体均值模型试验区域上的总体均值用试验点集P={x1,···,

xn}

上的响应值的样本均值估计Koksma-Hlawka

不等式：其中V(g)是函数的全变差，D*(P)为设计的星偏差。均匀设计可以满足要求，且稳健。75.3均匀性度量例试验区域:Cs=[0,1]s，P={x1,,xn}:Cs

上的n

个点.8D(XP)在X的行交换或列交换是不变的，即改变试验点的编号，或改变因素的编号，不影响D(XP)的值；若将XP

关于平面xj

=1/2反射，即将XP

的任一列(x1j,···,xnj)′变为(1−x1j,···,1−xnj)′，则它们有共同的D(XP)；D(XP)不仅能度量XP

的均匀性，而且也能度量XP投影至Rs

中任意子空间的均匀性。满足Koksma-Hlawka

不等式(5.8)；易于计算；与其它的试验设计准则有一定的联系，例如混杂、正交性，平衡性，等等；均匀性度量的要求：9Lp-星偏差其中，Fu(x)=x1…xs

为Cs

上的均匀分布函数，x=(x1,…,xs).

FP(x)表示设计P={x1,···,xn}的经验分布函数。10局部偏差函数其中则Lp-星偏差为11[0,x)

中的局部偏差函数xJx比率

=5/30 =0.167面积

=0.1615xJx比率

=9/30 =0.3

面积

=0.2875xJx比率

=6/30 =0.2

面积

=0.243812L2–星偏差Lp-

星偏差的缺点:

Dp(P)

旋转不可逆，原点

0

占有非常重要的位置没有考虑投影的均匀性13(a)D(P)=0.1611(b)D(P)=0.1500(c)D(P)=0.1411(d)D(P)=0.1389设计的旋转D(P)

14其中，需在

Rs

的所有低维上求和

u :{1,2,,s}

的非空子集;

xu :

X

在

Cu

上的投影；

Jx :预先由

x

确定的区域;

Jxu :Jx

在

Cu

上的投影。不同的

Jx

会导致不同的偏差定义.改进的偏差15中心化偏差对应的

JxxJx比率

=5/30 =0.167面积

=0.1615比率

=2/30 =0.0667面积=0.0478xJx比率

=2/30 =0.0667面积=0.0741Jxx16可卷偏差

(WD)17再生核希尔伯特空间希尔伯特空间：完备的内积空间核函数：18可分核：再生核19再生核希尔伯特空间20偏差定义：在偏差的定义中，目标函数取X上的均匀分布Fu。此时，偏差(5.25)式的具体计算公式为21(a)中心化偏差:22(b)可卷偏差:23(c)离散偏差:式中，当xij

=xkj

时δxijxkj

=1，否则为0。24(d)Lee偏差:式中βijk

=min{|xik−xjk|,1−|xik

−xjk|}25例5.2不同5水平设计的偏差的比较表5.2.两因素五水平的设计不同的偏差值2627例5.3不同6水平设计的偏差的比较表5.3.两因素六水平的设计28例5.3(续)295.4

均匀设计的构造均匀设计的基本要素因素个数和试验个数：s，n试验范围：超矩形、单纯性均匀性测度：某种偏差试验设计空间：全体试验设计PX

30均匀设计给定试验的诸元素(n,s,X)及均匀性测度D，若一个设计P∗∈PX

具有最小的偏差值，即D(P∗)=minD(P),则称P∗为在(n,s,X,D)下的均匀设计，或简称为均匀设计。求解均匀设计的方法分类：(A)理论求解；(B)优化数值求解；(C)求近似解。31均匀设计的理论解单因素试验中心化偏差可卷偏差不同的偏差，其相应的均匀设计可以不同也可以相同32均匀设计的近似解U

型设计

若

n×s矩阵U=(uij)中第j列的取值为1,2,···,qj

且这

qj

个数出现的次数相同，该设计称为U型设计，记为

U(n;q1,···,qs)。33缩小设计空间利用U型设计，寻求均匀设计的设计空间可大大缩小如下：Ũ(n;ns)为设计U(n;ns)的导出矩阵，即U

=

(uij)

=

U(n;ns)Ũ(n;ns)

=(xij)其中

xij

=(uij-0.5)/n34给定n和s，若设计的第一列选取为自然顺序1,2,···,n，则总共有(n!)s−1

个U型设计U(n;ns)，即使n和s不算太大，(n!)s−1

也是难以用穷举法来寻找均匀设计。好格子点法(glpm)

是有效的quasi-MonteCarlo方法

(N.M.Korobov(1959),E.Hlawka(1962),H.Niederreiter(1978),L.K.Hua,Y.Wang(1981),J.E.N.Shaw(1988),K.T.FangandY.Wang(1994).5.5好格子点法35Glpm

构造

Un(ns)

的步骤寻找备选的正整数集

Hn={h:h<n,gcd(n,h)=1.}，若Hn

中元素个数不小于s，转步骤2，否则glpm

无法构造需要的均匀设计，退出程序在

Hn

中选取

h1,h2,,

hs,产生一个

n

s

的矩阵

U=(uij)，

其中

uij=ihj

(modn),

其中模运算要求1

uij

n.记

U

为

U(n,h),其中

h=(h1,

,hs)称为

U

的生成向量.记

Un,s

为所有

U(n,h)

的设计.在Un,s

中选择使得偏差值最小的设计即为所求的(近似)均匀设计Un(ns).12336例

5.5对于

n=15，s=2

g.c.d{1,15}=1

g.c.d{2,15}=1

g.c.d{3,15}=3xxxxxxx37在中心化偏差意义下，当生成向量

h*=(1,11)时，U(15,h*)的偏差最小，故得均匀设计

U15(152).h*=(1,11)CL22=0.00160038好格子点法的注意事项生成向量中元素与n互质的的目的是使得向量hj

=(hj,2hj,···,nhj)(modn)构成{1,2,···,n}

的一个置换；给定n，其唯一的素数分解为，与其互质的正整数的个数为欧拉函数ϕ(n)欧拉函数ϕ(n)的个数多少对最后选择出来的设计的均匀性有很大影响。39例当

n=6,H6={1,5}，其元素个数只有2.则好格子点法得到的设计

U6(62)均匀性不好.因此,由

H6

更不可能得到U6(6s),s>240修正的好格子点法寻找备选的正整数集Hn+1在Hn+1

中选择s个不同元素

h=(h1,h2,,

hs),产生一个

(n+1)

s

矩阵

U(n+1,

h1,,

hs).

去掉第n+1行得到设计U(n,h),其全体记为

Un,s

.选择

使得

有最小偏差

则U*

为(近似)均匀设计.123U*U*

Un,s

41例(续)当

n=6,得

H7={1,2,3,4,5,6}修正的好格子点法U6*(62)好格子点法U6(62)所有

U6(6s),s

6的近似均匀设计可用H7

得到.近似均匀设计

U6(62)42方幂好格子点法构造

Un(ns)

的步骤寻找正整数集

An,s

={a:a<n,gcd(a,n)=1,且a,a2

···,as

在

模

n

的意义下互不相同}.对每个a∈

An,s，按下面方法生成

U型设计Ua

=(uija):

uija

=iaj−1(modn),i=1,···,n;j=1,···,s.这里同余运算是把数减去n的适当的倍数后落入区间[1,n]。在所有的Ua

中选择a∗

∈

An,s

使得Ua∗

的具有最小的偏差。则Ua∗

即为(近似)均匀设计Un(ns)。12343切割法构造

Un(ns)

的步骤初始设计：设Up

=(uij)为试验次数p

的均匀设计Up(ps)，其中p>n，最好p≫n，且p

或p+1为素数。行排列：对于

l∈{1,···,s}，矩阵

Up

的各行按照第

l列的元素的大小重新排序，使得第

l列变为1,···,p，变化后的矩阵记为Up(l)=(uij(l))。切割：对于m=1,···,p,设n×s

的矩阵Up(l,m)=(uij(l,m)),其中12344切割法构造

Un(ns)

的步骤新设计空间：对矩阵Up(l,m)

每一列的n个元素，根据其大小从最小到最大分别重记为1到n，重记后的矩阵为U(l,m)，则为U型设计。共有ps

个这样的U型设计。推荐的设计：对于给定的偏差准则

D，比较这

ps

个

U型设计U(l,m)，偏差最小的设计即为欲求的(近似)均匀设计

Un(ns)。45注：步骤1的初始设计可以用好格子点法构造，或当n,s较大时，用方幂glpm

构造。步骤5中的偏差可以为中心化偏差、可卷偏差，等等。45例5.6初始设计任选一列，如第三列，排序取m=13，应用(5.36)式得到构造一个均匀设计U9(93)46例5.6(续)按该列的大小关系，从最小到最大分别重记为1到9，所有得到的设计中，U(3,13)有最小的中心化偏差0.1027，比直接用glpm

得到的偏差0.1044要小。47切割法的优点给定一个初始设计Up(ps)，可以找到多个不同n和s的近似均匀设计Un(ns),n<p；用切割法构造的设计比直接用glpm

构造的设计均匀性更好；该设计不太依赖于偏差准则。48大试验次数的设计双因素试验

Fibonacci序列{Fn}：F0=0,F1=1,Fn+1=Fn

+Fn−1(n>1)。则对于试验次数为Fm(m≤3)的两因素试验，由(1,Fm−1)作为生成向量的设计具有很好的均匀性。多因素试验

设a1,···,as

为s个不同的正整数，记Fi

=ai(i

=1,···,s),Fm

=Fm−1+Fm−2+···+Fm−s+1(m>s)，则称序列{Fm}

为广义Fibonacci序列，则试验次数为Fm(m>3)的设计可由(1,Fm−1,Fm−2,···,Fm−s+1)作为生成向量构造，对较大的n，所获的设计具有较好的均匀性。49设

Un,s

为全体

U(n;ns)的设计.我们寻找

U*

Un,s

使其在给定的偏差准则D

下达到最小,即5.6随机优化法50优化问题的困难:

即使对于一般大小的n

和s，备选集是很大的有限集合。

对于该问题，传统的优化算法失效.局部极小值可能成百上千个.遗传算法门限接受法模拟退火因此，需要寻找更佳的优化算法

U*.51模拟退火算法的思想当新设计优于当前设计，即新设计的偏差低于当前设计，接受新设计，然后把新设计作为当前设计；若新设计比当前设计差，则用一定的小概率接受新设计门限接受法：当新设计比当前设计差的不太多时，直接接受新设计，换句话说，其判断准则采用硬门限，而不是用一定的概率。52设

f(x)为定义在

G

(

Rs)

上的函数.寻找

x*

使得对于任意x

G，其邻域

N(x)需事先给定。53门限接受法的示意图545.7均匀设计的应用确定因素：对于具体问题，选择因素和相应的试验范围、水平数；试验方案：选择相应的均匀设计安排试验方案；作试验：按试验方案进行试验，得到相应的响应值；建模：根据因素的取值及相应的响应值建立合适的模型；预报：根据建立的模型给出响应和因素间的关系，同时寻找其极值点；追加试验：具体步骤：555.7.1仅有定量变量的试验例5.7.在某化工的合成工艺中，为了提高产量，试验者选了四个因素：原料配比(x1)，反应温度(x2)，反应时间(x3)及加碱量(x4)，试验区域为[1,5.4]×[5,60]×[1,6.5]×[15,70].(5.49)每个因素均取了12个水平：原料配比(mol/mol):1.0,1.4,1.8,2.2,2.6,3.0,3.4,3.8,4.2,4.6,5.0,5.4;反应温度(℃):5,10,15,20,25,30,35,40,45,50,55,60;反应时间(h):1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5;加碱量(ml):15,20,25,30,35,40,45,50,55,60,65,70.全面试验次数：124=207365657多项式回归模型(逐步回归法)一次模型

ŷ=0.063+0.004x2+0.046x3,

R2=86.9%，s2=0.0023二次模型ŷ

=0.2259−0.0507x1+0.0096x12+0.0010x2x3

+0.00035x3x4

R2=98.5%，s2=0.0003，二次中心化模型没有非中心化的二次模型好58统计诊断59方差分析表表5.7.例5.6的方差分析表60寻求更佳的工艺条件不难求得，当ŷ在

x1=2.6406,x2=60,x3=6.5,x4=70,达到极大值，ŷ

=0.7082,它显著地高于12次试验的最高y值0.6096。应在最优点做追加试验，假如通过追加模型，说明该模型是合理时，由于x2,x3,x4

的取值都在其边界点，故需要扩大原来设定的取值范围，再做追加试验。61混合水平625.7.2含定性变量的试验伪变量方法设定性变量x的值来自于不同的类C1,···,Cq，则对x的数量化需要q−1个伪变量z1,···,zq−1：63回归模型广义线性模型二次模型变量筛选方法中，变量z1,···,zm

不