弱相关试验设计（修订）

贺深泽 (Email：[email protected] )

整理说明

本文原稿独立作于2005年，发表于《数学的实践与认识》（Vol.39，No. 3, 99-107）。 经查，零相关与弱相关试验设计的概念与国际上一些数学家的 orthogonal Latin hypercube 及 nearly orthogonal Latin hypercube 类似。和文献相比，我的表述不大相同。例如，正交性度量引入临界值。由于特殊原因，发表时未修未校。现予修订，整理成 HTML 版本，予以公布，也为方便本博客读者。

HTML 写数学表达式很不方便。有兴趣的读者请查找原文。

摘 要

本文把选择交换向量分量的方法引入蒙特卡诺法，搜索一个与已知向量零相关的向量的速度大大增加， 构造零相关阵列成为可能。 运行数 n 在范围 4≤n≤33 ( n≠4k+2)内的 n×(n-1) 维弱相关阵列已经被构造出来，其中包含一个零相关子阵。 它们的试验点均衡分散，回归系数估计简单且相关性小。 文中给出了一个包含 7×3 维零相关子阵的7×6维弱相关矩阵及其正交化矩阵的例子。

关键词

试验设计，正交设计，零相关设计，弱相关设计，选择向量分量交换法。

Abstract

After the method of choosing to exchange vector components is introduced into the Monte Carlo method, the speed of searching a vector which is zero-correlation with a existing vector is greatly increased, thus the construction of zero correlation arrays is become a possible. The n×(n-1) weak correlation arrays have been constructed, for the runs n at the range 4≤n≤33( n≠4k+2), a zero correlation sub-array be included in them. Their experimental points are well-balanced and the estimation of the regression coefficients are simple and the correlation between the estimates is the weak. An example of 7×6 weak-correlation array which included a 7×3 zero-correlation sub-array is provided in the paper.

Keywords:

Experimental designs, Orthogonal designs, Zero-correlation designs, Method of exchanging Selectively components of a vector.

一. 引 言

设一个 n×p 维设计矩阵为 X=(x₁,…,x_m), 其中 x_i=(x_1i,…,x_ni)^T， 统计学定义其相关矩阵为 R=(r_ij),其中，

r_ij =cov(x_i,x_j)/(σ_iσ_j) (i,j=1,2,…,m) ------ (1.1) (注：受语言限制，本文中变量的均值用变量右上角加“-”号表示。类似，回归系数估计值右上角加“^”号表示。） x_i与x_j的协方差cov(x_i,x_j)由 ∑ⁿ_k=1(x_ki-x^-_i)(x_kj-x^-_j)/(n-1) 估计， x_i 的方差 σ²_i由&Sum;ⁿ_k=1(x_ki-x^-_i)²/(n-1)来估计。 两个n 维向量的相关系数 r 小于相关系数临界值 r_α(n-2)，则说它们的相关性置信水平不小于α (整理时注：即其置信概率大于 1-α, 统计量 p=1-α 常常被称为 P 值)。 相反，如果r>r_α(n-2) 则它们的相关性置信水平达到或超过了α。 α 没有小到一定的数量不能认为相关性认定有意义。如果相关矩阵 R 中绝对值最大的非主对角线元素 mcc=max{ |r_ij| , i≠j, i,j=1,2,…,m} 大于相关系数临界值r_α(n-2)^[5,6]，则该设计至少有两个列向量相关性的置信水平小于α， 此时称该设计的相关性置信水平小于α。 通常，如果 α>0.25 (即，p≤0.75)，便认为该设计为弱相关设计。 同理，如果 α<0.25 (即，p≥0.75)，则说该设计为强相关设计，α=1(p=0) 时为零相关设计。 在具体工程中 P 值的水平根据工程的要求确定。

暂时给设计X添加一个1 向量 x₀=1 做第一列，仍记为 X。统计学^[2,3,5,7] 有如下的结果，对于线性模型

y = b₀+ b₁x₁ + b₂x₂+...+b_m x_m ------(1.2) 最小二乘系数估计 β^{^} 由正规方程（矩阵式） Ab=B ------(1.3) 的解给出。这里，A= X^TX，B= X^TY。 只要行列式 det A≠0，就能得到唯一的估计 β^{^}=A^-1B ------(1.4) 而且，估计β^{^}的协方差矩阵为 Var(β^{^})=A^-1σ². ------(1.5) 此处的σ²由(σ^)²=∑ⁿ_i=1(y_i-y^_i)²/(n-m-1) 来估计。 这里，y^_i 记 y 在第 _i 个试验点的预报值。 这意味着，一般来说，回归系数估计值之间有相关性，除非 A^-1 是一个对角矩阵。A^-1 的非对角元素如果很微小，则相关性微弱。

当X是正交的，即, 列向量的内积 (x_i,x_j)=0 (i≠j；i,j=1,2,…,m)，且 x^-_i=0, x^-_j=0, 则，A 是对角阵。

对原始设计的变量x_i做标准化线性变换

x_ki=(x_ki-x^-_i)/σ_i (k=1,2,…,n; i=1,2,…,m)------(1.6) 将移动试验区域到以原点为中心的位置并调整比例,不改变点与点之间的位置关系, A将具有形式 其中A₁ 是 n 维行向量。每个设计都对应一个相关矩阵。只要设计的变量之间有相关性，参数估计之间就会有相关性。变量之间的相关性越小，参数估计之间的相关性越小。如果变量之间没有相关性，则参数估计之间没有相关性。逐步回归分析^[7]允许变量之间有一定的相关性，有能力从回归变量的集合中选取显著变量子集。通过设置蜕化允许值界限来防止矩阵蜕化，并阻止高度相关的两个因子同时被选取。倘若某两个变量完全相关，A 不满秩，(1.4) 不能被唯一地确定。某些方法，例如岭回归或直接从叉乘和矩阵获取，可以得到参数估计。两个相关的变量中编码在前的变量独占了两个变量的方差被保留，后一个弃用。对样本的拟合看起来天衣无缝。但可质疑，被保留的变量一定是显著的吗，被弃去的一个一定是不显著的吗？它对变量的功能与效应的解释常常是不正确的，因而在整体上，它的预报偏差很大。如果某个变量与其他变量之间的相关性很强，同样会面临矩阵蜕化的局面，回归系数估计不稳定。

如果X是正交的，R将是单位矩阵，A₁为0 向量。如果 X 的任何两个变量之间零相关，则 R 也会是单位矩阵，A₁也是0 向量。这就显示了零相关设计的重要意义。

以 (1,2,…,n)作各个变量的基本水平设计，对 n>2，田口玄一以n²个试验构造出了n²×(n+1)维设计矩阵，称为正交表^[8,6]。 按代数学观点，当应用于定量因子时，它们是零相关的，并不正交。应用这类矩阵作连续量的回归设计时，当以原点为中心关于原点对称地进行水平设计时，任何两个列向量的内积(x_i,x_j)=0, i≠j; i,j=1,2,…,n+1，满足代数学上的正交性定义。这类表的最大特点是具有均衡分散齐整可比的特点，每个列向量的每个水平的重复次数相等；任意两个列向量a,b的水平交互配合所构成的数对(a_i,b_j) 具有相同的次数。这是正交设计的第二种定义方式。用于非数字和不连续量的试验，不但能够估计出每个列的效应，还可以估计出每个列的每个水平的效应。

正交表也存在一些问题，例如用作连续量的回归设计时，二水平试验的试验点密度太低，单由一组试验所能得到的结果很有限。初学者不会处理，往往把过程一律当作线性过程处理，试验次数虽然少，分析推断得不到正确的结论，很容易陷入困局。想“撒大网，捉大鱼”，结果没有捉到鱼。多水平试验试验次数太多，难以接受。

方开泰的均匀设计试图解决正交设计的这一问题 ^[9]。现有两类均匀设计表，由方开泰在70年代末用同余乘法计算出来的设计表(以下简称同余表)和由马长兴用蒙特卡诺法计算出来于1999年发表的设计表^[10](以下简称新均匀表)。同余表的偶数表与奇数表只有最后一个试验不一样，由奇数表删除最后一个试验即得到对应的偶数表。同余表属于高相关性的设计。偶数表的秩不超过列数的一半，它们最多只能安排列数的一半那么多的试验因子。多则必包含完全相关的列向量。奇数表的高相关列数也大于列数的一半。矩阵规模越大，相关性越高。新均匀表在相关性方面有了很大的改进，但其结构不尽合理，某些设计表中还存在高相关的情况。

本文首先讨论在田口玄一正交表之外零相关设计的构造方法及其基本性质，然后构造弱相关矩阵。

二．零相关矩阵存在的条件与构造方法

2.1 零相关矩阵的存在条件

用h=(1,2,…,n)^T记基本列向量。h的分量的不同排列方式有n！种，都是n维实欧氏空间Ε_n中的向量，且都是Ε_n的格点。它们构成一个向量集合S。

设一个零相关矩阵由 p 个列向量构成 x_i=(x_1i,…,x_ni)^T∈S,(i=1,2,…,m),x₁=h。r_ij= 0. r_ij=0 等效于叉乘和等于0,

∑ⁿ_k=1(x_ki-x^-_i)(x_kj-x^-_j)=(x_i,x_j)-cf =0, (i≠j；i,j=1,2,…,m) ------(2.1) 这里cf=nx^-_ix^-_j。这个方程在S中不一定有解。我们有结果，

定理：零相关矩阵存在的必要条件是 n 为大于 3 的奇数或 4 的倍数。当 n=4k+2（k为非负整数）时不存在零相关矩阵。

证明：要想(2.1)在S中有解, cf 必须是整数。因为诸 x 的分量都是非负整数，它们的内积一定是非负整数，如果cf 不是整数，(2.1) 恒不能成立。因此，当cf 不是整数时，零相关矩阵对任何 p (列数) 不存在，此时(2.1)在S中无解。

由于cf= nx^-_ix^-_j= n(1+ n) ²/4，分母为4， 若 n 为 4 的倍数，4 是 cf 的分子的一个因数，cf 必为整数。 若 n 为奇数，(n+1) 必为偶数，(1+n) ²必包含因子4，cf 的分子中必包含因子4，分母是分子的约数。

n 不能是 3，当 n=3 时，S 中只有 6 个向量，可以直接验证，任何两个的相关系数不等于零。

综上所述，零相关矩阵存在的必要条件是 n 为大于 3 的奇数或 4 的倍数。

当 n 为偶数但不是 4 的倍数时必为 2 的奇数倍数，(n+1) 必为奇数，平方后仍为奇数。此时，cf 的分子必不是 4 的整数倍数，因而必不是整数。因此，此时零相关矩阵不存在。

这个条件不充分。例如，当 n=3 时，cf(=12) 是正整数，零相关矩阵不存在。对于3< n <34, n≠4k+2, 零相关矩阵都已经被构造出来。但尚不能证明对于更大的 n 一定存在。

2.2 零相关矩阵存在的不唯一性

在 S 中与h零相关的向量应当满足条件

(h,x)- cf = 0 ------(2.2) 这里有 n 个未知数，只有一个方程，有无穷多解。指定n-1个变量的值，所得到的解不一定在 S 中。 因而不能用解不定方程的方法来求解。然而，容易验证，当 n=5,7,8 时，与h 零相关的向量分别有6个，184个和936 个。因此，n×2维零相关矩阵如果存在，不是唯一的。第三节中将进一步阐明零相关矩阵的多态性。

2.3 零相关矩阵的序贯构造方法

即使零相关矩阵存在，矩阵有多大是未知的。因此，必须序贯地逐列构造。 假定已有 p 个列为c_i=(c_li,…,c_ni)^T,(i=1,2,…,m)， 其中 c₁= h。待求新列 x =(x₁,…,x_n)^T应当满足的条件是

(c_i,x)-cf = 0, ( i=1,2,…,m) ------ (2.3) 用S₁记S 中与h(即c₁)为零相关的向量的子集(即(2.2)解的集合)，x∈S₁。S₁的成员相互之间的相关系数不都是 0，因而不能用简单的方法从S₁ 得到新列向量。设S_i 是与c₁,…,c_i同时零相关的向量的集合，x∈S_i。i 越大，交集越小。最大化的零相关矩阵是S 的一个能使全部列向量相互零相关的最大子集。下面给出三种序贯构造方法。

2.3.1 全排列筛选法

1．c₁= h；(整理时注：新版本(2017)的第一列均采用 h 的一个随机排列。)

2．对S的成员逐一检验，如果x∈S满足(2.2)，以x为第2列；

3．假设已找到 i 个列，如果S中的一个成员满足(2.3)，则以该列为第 i+1列。依此类推，直到不存在一个向量满足(2.3) 时，结束选取。此时称它已经零相关饱和。

随着 n 增大，计算工作量爆炸性地增加。全排列筛选法所需时间太长。

2.3.2 随机法

随机法通常叫蒙特卡诺法。随机地从 S 中选取一个向量 x，如果它与各已知列向量皆为零相关便以它作为所求的一个新列向量。 否则，随机地从 S 中另取一个向量，重复上述过程。

即使是 m=2，任选一个向量满足(2.2)其概率也非常小。当 n=7 时,一次选中的概率为184/7！≈0.0365；n=8 时,一次选取的概率大约为936/8!≈0.0232。更不用说满足(2.3)了。

2.3.3 选择向量分量交换法

引理: 设 c=(c₁,…,c_n)T∈S, x=(x₁，…,x_n)^T∈S , 记 g=(c, x)。交换 x 中的两个分量，内积的增量为 c 与 x 对应分量的差之积，两对分量的差同号时，符号为负，异号时符号为正。

设 x 中被交换的分量的下标分别为 i,j(不失一般性，设 i<j),x' 为交换得到的新向量，新老内积之间的关系为 g’= g – c_ix_i-c_jx_j+c_ix_j+c_jx_i = g – (c_j–c_i)(x_j-x_i) ------(2.4) 这意味着当 c, x 中的两对分量的差同号时，交换后内积减小，异号时内积增加。

由于交换 x 中任两个分量之后，a=|(c, x)-cf|或者减小，或者增大。因而可以通过选择分量进行交换，使 a 减小。 在有限次交换之后，a 可以持续减小，甚至收敛于0。随机地选择一个向量恰好 a=0 的概率非常小，经过若干次选择分量交换使之等于 0 的概率就大得多。

假设已找到 k 个列向量，第 k+1 列向量与其前 k 个列向量有 k 个 a_i=|(c_i, x)-cf|, 交换 x 中两个分量后，某些a_i可能减小，某些a_i 可能增大。选择使 a_i 一致地变小收敛于 0 的分量进行交换，速度更快。

三．零相关矩阵的基本性质

3.1可堆叠

维数相同的两个零相关矩阵上下堆叠之后仍然是一个零相关矩阵。

设 A=(a₁,…,a_p)和 B=(b₁,…,b_p)是两个n×m 维零相关矩阵， 堆叠构成一个2n×m 维矩阵C=(c₁,…,c_m)。容易得到

(c_i,c_j) – cf_c = (a_i,a_j)- cf_a+ (b_i,b_j)- cf_b ------ (3.1)

维数相同的不同零相关矩阵堆叠后，试验点密度增加，均匀性互补。

（整理时注：这个命题不够准确，详见 堆叠法 )

3.2 不可拼接

饱和的零相关矩阵拼接（或嵌入）任何一个属于S的列向量都不可能再是零相关矩阵。如果存在与它的已有列向量零相关的列向量，则它就不是零相关饱和的。因此，饱和的零相关矩阵不能拼接使用。

3.3 剪裁性质

一个零相关矩阵的部分列向量所构成的矩阵仍然是一个零相关矩阵。但部分行向量所构成的矩阵不是零相关矩阵。

3.4 多态性

行(列)可移动 零相关矩阵的行（或列）向量平行移动、交换之后仍然是零相关矩阵;

列的多选性 在零相关饱和之前，列有多选性，选择不同，后面的列向量结构不同。甚至矩阵的列向量数目不同。例如，n=8 时可能只有 3 个列向量，也可以有 4 个列向量，n=9 时可能找到 4 个列向量，也可能是 5 个列向量, 取决于搜索方法、搜索深度和选取路径。

多种组合方式 利用列向量的移动、矩阵的堆叠和剪裁等性质可以组合得到多种形态的零相关矩阵。

3.5 正交性

当一个零相关矩阵所有变量的水平设计等距离且均值x^-_i=0,(i=1,2,…,m)，此时设计是正交的。 如果水平设计为等距离的而均值不等于0，则设计不是正交的, 但线性变换z_ij=x_ij-x^-_i之后为正交。 如果变量的水平设计不满足上述条件，有可能破坏零相关特性，更不满足正交性条件。

3.6 均匀性

在概率论中^[2,4]有结果：n个随机变量ξ_i(i=1,2,…,n)在区域 D={(a_i,b_i),i=1,2,…,n}上的分布如果是均匀的， 每个随机变量ξ_i 的数学期望Eξ_i=(a_i + b_i)/ 2，方差是一个 n 阶对角矩阵，对角线上的元素依次为各随机变量的方差。 对于离散随机变量，这个命题其逆不真。但如果是随机取样，在x^-_i=(a_i + b_i)/2,(i=1,2,…,m) 的条件下，零相关设计不均匀的概率非常小。在这一意义上，零相关设计是均匀的。（整理时注：考虑到均匀的统计学特性，使用术语均衡性更恰当。

零相关矩阵正交化^[3]后，列向量的长度不变，这就意味着试验区域的几何形状和体积不变，试验点充满了整个试验空间。因此，在欧氏空间的几何分布的意义上，它们也是均衡的。

四． 弱相关设计矩阵的构造方法

有两种情况引起构造弱相关设计的要求：当 n=4k+2 时不存在零相关矩阵，可以用弱相关矩阵来近似；当零相关列数太少，可以补充一些弱相关列来满足设计的需要。

弱相关矩阵可以看作为含有微弱误差的零相关矩阵，微弱的相关性可以忽略，全部回归系数的最小二乘估计相关性微弱。

设一个弱相关矩阵的 p 个列向量为c_i=(c_1i,…,c_ni)^T,(i=1,2,…,p),c₁= h，若x=(x₁,…,x_n)^T∈S，与已有 p 个列向量的 p 个叉乘和的绝对值分别为

r_k(x)=|∑ⁿ_k=1c_ikx_i - cf|, k=1,2,…,m------(4.1) 求S 中(4.1)的极小解 min max{r_j(x),j=1,2,3,…,k } ------(4.2) 即求使全体r_k(x)小于一个足够小的数值的向量 x。 同零相关矩阵一样，可以序贯地构造。 如果存在零相关矩阵，在零相关矩阵的基础上，继续采用前面的方法，逐步追加新的列向量； 如果不存在零相关矩阵，从第 2 个列向量开始追加新列。可以得到n×(n-1)维弱相关矩阵。

本文方法构造的弱相关矩阵各级子阵的最大相关系数(mcc)逐列递增不减。要使设计获得小的相关性，可以空闲出右边的某些列，直至只用零相关列。

五．弱相关矩阵示例

以n=7为例，弱相关矩阵 W7 (其中包含一个7×3 的零相关子阵) 示于表 1，其相关矩阵示于表 2。

表1. W7 表2. W7 的相关矩阵(上半部分)

有时希望构造 n×n 弱相关方阵，正交化后一定是相关的。

六．弱相关设计的分析计算方法

有条件时，用回归分析程序处理数据，可以参考各种有关回归分析的专著 ^[3,5,7]。

6.1 直接比较

由于设计是均衡的，把从 n 个试验的结果的直接比较中得到的较好的条件认定作为总体的一个优化条件是有优势的。 它不是推理的结果，具有现实性，靠近总体的最优值的可能性很大。在这一意义上，弱相关设计用于含有非数字和不连续量的试验也是可行的。

6.2 线性回归分析

由于变量间弱相关，回归系数估计值之间弱相关，从叉乘和矩阵即可以得到回归系数的估计值，微弱的相关性可以忽略。 回归计算时，微弱的相关性会完全消除。微小的回归系数删除之后，不必重新建模，重新计算。因此不必做逐步回归。 如果事前的设计为正交的，逐步回归分析结果与从叉乘和矩阵中得到的结果相同，还提供了回归常数的准确估计。

为演示计算过程，下面模拟一个线性过程。假设

y = 71.8 + 5.12x₀ – 0.73x₁ + 9.62x₂ –0.0083x₃ -----(6.1) 使用w7，用(6.1)作试验数据发生器。表3 为其模拟样本的叉乘和矩阵，最底下的一行紧凑地安排了回归系数估计值。这里的w7与上节中的w7稍有不同，正好展示了w7的多态性。 表3. w7模拟数据的叉乘和矩阵 由于叉乘和矩阵中与试验变量有关的非对角元素大多为0或接近于0，可以按表 4 的方式估计回归系数，要计算的只有那几个粗体字。 如表 4 所示，由于设计方案相关性非常微弱，接近正交，均值也提供了回归常数的准确估计。 由于x₁与x₃有微弱的相关性, 这两个变量的回归系数有微弱的相关性, 估计值与原值有微小的误差。 可以相信，如果采用弱相关矩阵w8(它有四个零相关列)，结果会更好。因此，对于n=4k+2 的那些弱相关设计，宁可用比它大 1 或小 1 的那个设计。

表4. W7正交点上模拟数据的分析计算表

6.3 非线性回归分析

如果过程为非线性的，为了建立非线性回归模型，需要导出新变量。 导出变量通常是原试验变量的一个函数，它与原变量必然会存在一定的相关性。 应该用相关矩阵检查回归模型，如果包含相关或甚高相关的回归变量应该调整， 删除相关性太高的回归变量让位给另一些回归变量，以利提高效率。 非线性回归一定要用向后逐步回归计算，否则，不能有效地消除回归系数之间的相关性。分析推断的效果不佳。

工艺优化最适宜使用梯度法，在小扰动条件下，线性模型回归系数提供了梯度方向的良好估计。

七． 结论与存在的问题

为了说明强相关设计存在的问题，可以用如下的模拟。表5 是一个强相关设计方案（简称Q），这个设计的模板是均匀表 u*6(6⁴)。 试验数据发生器仍为(6.1)。表 6 为试验样本的相关矩阵，表 7 为其叉乘和矩阵。因为矩阵不满秩，回归扫描出错。 表 7 表明，由叉乘和矩阵得到的四个效应的估计值(最底下一行)与真值相差很大，完全相关的x₀与x₃的效应估计绝对值相等，符号相反。 但容易验证，以下两个方程

y = 71.8 + 5.12x₀ – 0.73x₁ + 9.62x₂ ------(7.1)

y = 71.8 – 0.73x₁ + 9.62x₂ – 5.12xx₃------(7.2) 对表 5 中试验数据的拟合都没有偏差。但用于表 4 中试验数据时，偏差则大得不能接受。说明在高相关设计中不能仅根据对样本的拟合好坏来决定接受或拒绝一个拟合函数。在几何上，表 5 的试验点分布在一条“直线”上，不包含“线”外的信息。这也证明了它的极端的不均匀性。 表5. Q方案及模拟试验数据 表6. Q样本的相关矩阵(上半部分) 表7. Q样本的叉乘和矩阵(上半部分)

目前已经在单台微机上构造出了直到 33×32 维弱相关设计（整理注：2017 版已经更新到 n=37,而且，正交列数普遍增加）。 试验次数少，效率高。模拟展现了其易使用，计算简单，估计误差小的特点。预示其有着良好的应用前景。

现有结果，当 n>9, 零相关子阵没有超过 5 个列向量。究竟应该含有多少个零相关列向量，需要继续研究。

参 考

1．George E.P.Box, J.Stuart Hunter & William G.Hunter Statistics for Experimenters (Second Edition ) Wiley-Interscience (2004)

2．J.S.Milton, Jesse C.Arnold Introduction To Probability and Statistics (Third Edition) Irwin McGraw-Hill （1995）

3．T.Hastie , R.Tibshirani & J.Friedman The Elements of Statistical learning （Springer 2001）

4．复旦大学数学系主编 《概率论与数理统计》（第二版） 上海科学技术出版社 （1961）

5．茆诗松 丁元 周纪芗 吕乃刚编著《回归分析及其试验设计》华东师范大学出版社（1981）

6．中国科学院数学研究所统计组《常用数理统计表》科学出版社（1974）

7． R.I.Jennrich，杨自强译 逐步回归《数字计算机上用的数学方法》卷Ⅲ(4)，科学出版社（1981）

8． 田口玄一 《实验计划法》第三版 丸善株式会社 东京 （1976）

9． 方开泰 均匀设计 应用数学学报 卷3（4）（1980）

10．方开泰，马长兴 著 《正交与均匀试验设计》 科学出版社 (2001)