CN105068972A - 一种拉丁超立方采样的扩展方法 - Google Patents

Abstract

本发明公开了一种拉丁超立方采样的扩展方法。假设已有样本大小为n的LHS，扩展方法的目的是得到样本大小为n+m的LHS新采样，同时新采样中最多地包含已有采样点。本发明的方法通过将LHS采样点之间的关系转化为简单无向无环图，得到邻接矩阵，再求解最大独立集来最多地保留已有采样点。实现步骤包括计算已有采样点在新采样结构中的分布，计算需要保留的采样点矩阵，计算新增采样点矩阵，最后计算扩展后的拉丁超立方采样矩阵。本发明能够尽量少的删除已有采样点，节约采样成本；同时生成新采样点，使得新采样点和保留的采样点所构成的采样样本，仍满足LHS结构。

Description

一种拉丁超立方采样的扩展方法

技术领域

本发明涉及数据信息处理领域，具体而言，涉及一种拉丁超立方采样的扩展方法。

背景技术

拉丁超立方采样(LatinHypercubeSampling，LHS)是一种全空间填充且非重叠的随机采样方法。全空间填充采样确保即使在没有详细的源函数特性的情况下，也可以得到该函数在整个设计空间的信息；非重叠采样则确保没有重复和多余的采样点。LHS也是一种随机采样方法，但是与随机采样不同的是，它产生的采样点在全局内是均匀的。传统拉丁超立方采样算法基于固定采样数，在应用中需要根据预先给定的采样数构造固定的采样结构和采样点，然后计算得到结果样本；若计算结果不理想则需要增加采样数，必须根据增加后的采样数重新构造新的采样结构并计算新的结果样本，无法利用原有的结果样本，降低了再次采样的效率。因此在已有采样结构的基础上扩展采样数，充分利用原有的计算结果样本是十分重要的。

拉丁超立方采样扩展(ExtensionofLatinHypercubeSampling,ELHS)问题是：已有样本大小为n的LHS，希望得到样本大小为n+m且满足LHS结构的采样样本。ELHS的目的是新的采样点中最大限度地利用已有采样点，这样可以在保证采样结构的基础上减少重新采样的时间和成本。

传统的拉丁超立方采样扩展方法首先构造样本大小n+m的LHS结构，并判断已有采样点在该结构中的位置，然后根据空白层的位置和个数产生新的采样，进而构成新采样。传统扩展方法中，新采样并不能保证是LHS，这是因为样本大小的变化导致LHS采样空间结构的变化，已有采样点可能不满足LHS结构，如图2所示，已有采样点在新的结构中有部分点落在同一区间内，不再是严格意义上的拉丁超立方采样。

发明内容

本发明提供一种拉丁超立方采样的扩展方法，针对LHS的扩展问题，在不破坏拉丁超立方采样结构的基础上最大限度地利用已有采样点，减少了重新采样的时间和成本。

为达到上述目的，本发明提供了一种拉丁超立方采样的扩展方法，包括以下步骤：

步骤1：计算已有采样点在新采样结构中的分布，包括计算原始排序矩阵映射到新采样结构中的排序矩阵原始随机数矩阵映射到新采样结构中的随机数矩阵其中排序矩阵的矩阵元素其中i＝1,2,…,n，j＝1,2,…,r，为向下取整；随机数矩阵的矩阵元素其中i＝1,2,…,n，j＝1,2,…,r，n、r均为自然数；

步骤2：计算需要保留的采样点矩阵首先根据构建采样点间的邻接矩阵A_n×n，将采样点间的关系用矩阵中的元素表示；然后对采样点1,2,…,n的每个组合进行全排列，全排列矩阵为B＝[B₁B₂…B_q]^T；取p＝1，根据B_P中排列组合的序号剔除对应采样点，即A_n×n中对应的行和列，得到新的邻接矩阵A'，如果A'不为零阵，则p＝p+1，重新删除采样点构建；如果A'为零矩阵，记B'_p为B_P中非零元素构成的向量，则B'_p中的元素就是剔除采样的序号，A'的维数记为n'，即保留n'个采样点；最后按照向量B'_p元素的值剔除和中对应的行，分别得到和其中邻接矩阵A_n×n的矩阵元素为

步骤3：计算新增采样点矩阵首先生成(n+m)×r维的新排序矩阵L'_(n+m)×r，将L'_(n+m)×r的每一列与相应列中相同的元素删除，将剩余元素构成的矩阵记为然后生成随机数矩阵由此计算累积概率值矩阵最后求得新增采样值矩阵

步骤4：计算ELHS采样矩阵根据步骤1、2、3得到新的采样值矩阵即为样本大小为(n+m)的LHS，其中采样值矩阵 $X_{(n+m)\×r}^{\left(new\right)}=\left[\begin{array}{c}{X}_{{n_{\′}}^{\′}\×r}^{\left(2\right)}\\ X_{(n+m-n^{\′})\×r}^{\left(3\right)}\end{array}\right],$ 其中分别由步骤2、步骤3求得。

进一步地，步骤2中所述的全排列矩阵为B＝[B₁B₂…B_q]^T，其中B的每一行都是1,2,…,n的一个排列，不足n个元素的最后用0补齐，且B_i中0元素的个数不少于B_j，1≤i≤j≤n。

进一步地，步骤3中所述的新排序矩阵L'_(n+m)×r的每一列都是随机生成的整数1到n+m的一个排列，随机数矩阵的矩阵元素r_i,j为[01]上均匀分布的随机数，累积概率值矩阵的矩阵元素的计算公式为

$u_{i,j}^{\left(2\right)}=\frac{l_{i,j}^{\left(2\right)}-1+r_{i,j}^{\left(2\right)}}{(n+m)}.$

进一步地，步骤3中所述的新增采样点矩阵的计算公式为

$x_i=F_i^{-1}\left({\mathrm{Prob}}_i\right).$

本发明提供的生成新采样点的方案使得新生成的采样点和保留的已有采样点仍能构成一个LHS；针对LHS的扩展问题，本发明在不破坏拉丁超立方采样结构的基础上最大限度地利用已有采样点，减少了重新采样的时间和成本；本发明对已有采样点和新采样点的个数没有要求，普遍适用于一般增量扩展、整数倍扩展等场合。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为变量维数为2，样本大小为5的拉丁超立方采样的示意图；

图2为拉丁超立方扩展导致的已有采样点在新采样结构中的分布的示意图；

图3为本发明将最大保留问题转化为图论中最大独立集问题所构造的简单无向无环图的示意图；

图4为本发明方法的算法流程示意图；

图5为依照本发明实施例的采用传统拉丁超立方采样扩展方法结果的示意图；

图6为依照本发明实施例的采用本发明方法扩展结果的示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的拉丁超立方采样的扩展方法将已有采样的最大保留问题转化为图论问题，将采样点映射为图中的顶点，采样点序号对应顶点序号，如果两个采样点中至少存在某个维度的变量落在同一区间内，则在图中将对应的两个顶点连接，即图中的一条边，构成一个简单无环图，如图3所示。这样已有采样点的最大保留问题就转化为图论中最少删除顶点使得图中边的数量为零，即图论中的最大独立集问题。本发明方法的数学描述及实现所采取的技术方案如下：

假设已有样本大小为n的拉丁超立方采样，样本矩阵为排序矩阵为累积概率值矩阵为随机数矩阵为需要删除n-n′个已有采样点，并生成n-n′+m个新采样点使得

$X_{(n+m)\×r}^{\left(new\right)}=\left[\begin{array}{c}{X}_{n^{\′}\×r}^{\left(old\right)}\\ X_{(n+m-n^{\′})\×r}^{\left(2\right)}\end{array}\right]$

构成一个样本大小为n+m的LHS。具体包括以下步骤：

步骤1：计算已有采样点在新采样点结构中的分布。

计算映射到新采样结构中的排序矩阵矩阵元素其中i＝1,2,…,n，j＝1,2,…,r，为向下取整；

计算映射到新采样结构中的随机数矩阵矩阵元素 $r_{i,j}^{\left(1\right)}=u_{i,j}^{\left(old\right)}\·(n+m)-(l_{i,j}^{\left(1\right)}-1),$ 其中i＝1,2,…,n，j＝1,2,…,r。

步骤2：计算需要保留的采样点矩阵

根据构建采样点间的邻接矩阵A_n×n，将采样点间的关系用矩阵中的元素表示为

对1,2,…,n的每个组合进行全排列，全排列矩阵为B＝[B₁B₂…B_q]^T，其中B的每一行都是1,2,…,n的一个排列，不足n个元素的最后用0补齐，且B_i中0元素的个数不少于B_j，1≤i≤j≤n；取p＝1，根据B_p中排列组合的序号剔除对应采样点，即A_n×n中对应的行和列，得到新的邻接矩阵A'；如果A'不为零阵，则p＝p+1，重新删除采样点构建；如果A'为零矩阵，记B'_p为B_p中非零元素构成的向量，B'_p中的元素就是剔除采样的序号，A'的维数记为n'，即保留n'个采样点；

按照向量B'_p元素的值剔除和中对应的行，分别得到和

步骤3：计算新增采样点矩阵

生成(n+m)×r维的新排序矩阵L'_(n+m)×r，其中L'_(n+m)×r的每一列都是随机生成的整数1到n+m的一个排列；

将L'_(n+m)×r每一列与相应列中相同的元素删除，剩余元素构成的矩阵记为

生成随机数矩阵矩阵元素r_i,j为[01]上均匀分布的随机数；

计算累积概率值矩阵矩阵元素

由公式求得采样值矩阵

步骤4：计算ELHS采样矩阵根据步骤1、2得到新排序矩阵新随机数矩阵新采样的累积概率值矩阵和新的采样值矩阵

$L_{(n+m)\×r}^{\left(new\right)}=\left[\begin{array}{c}{L}_{{n_{\′}}^{\′}\×r}^{\left(2\right)}\\ L_{(n+m-n^{\′})\×r}^{\left(3\right)}\end{array}\right]R_{(n+m)\×r}^{\left(new\right)}=\left[\begin{array}{c}{R}_{{n_{\′}}^{\′}\×r}^{\left(2\right)}\\ R_{(n+m-n^{\′})\×r}^{\left(3\right)}\end{array}\right]U_{(n+m)\×r}^{\left(new\right)}=\left[\begin{array}{c}{U}_{{n_{\′}}^{\′}\×r}^{\left(2\right)}\\ U_{(n+m-n^{\′})\×r}^{\left(3\right)}\end{array}\right]X_{(n+m)\×r}^{\left(new\right)}=\left[\begin{array}{c}{X}_{{n_{\′}}^{\′}\×r}^{\left(2\right)}\\ X_{(n+m-n^{\′})\×r}^{\left(3\right)}\end{array}\right]$

由此得到的即为样本大小为(n+m)的LHS。

为评估本发明方法的性能及应用，验证方法的有效性，并使本发明的目的、技术方案和有益效果更加清楚明白，下面结合具体实施例，并参照附图，对本发明进一步详细说明。

原始的拉丁超立方采样以附图1为例，即采样点维数为2，其中x₁服从[0,10]上的均匀分布，x₂是低限为5，众数为7.5，上限为10的三角分布，采样点个数为5。LHS的基本步骤为首先在每个维度的取值范围内等概率地分成5个相邻且不重叠的区间，每个维度等概率地随机选择一个区间，得到每个维度的分布函数F(·)为

F₁(x₁)＝0.1×x₁,0≤x₁≤10

$F_2\left(x_2\right)=\left\{\begin{array}{cc}\frac{{(x_2-5)}^2}{12.5}& 5\&le;x_2\&le;7.5\\ 1-\frac{{(10-x_2)}^2}{12.5}& 7.5<x_2\&le;10\end{array}\right.$

然后根据每个维度的分布函数F(·)的反函数求得采样值求得

x₁＝10×F₁(x₁)， $x_2=\left\{\begin{array}{cc}5+\sqrt{12.5\&times;F_2\left(x_2\right)}& 0\&le;F_2\left(x_2\right)\&le;0.5\\ 10-\sqrt{12.5\&times;(1-F_2(x_2))}& 0.5<F_2\left(x_2\right)\&le;1\end{array}\right.$

从而得到一组采样X＝[x₁x₂x₃x₄x₅]，其采样矩阵和相关矩阵为

$X_{5\&times;2}^{\left(old\right)}=\left[\begin{array}{cc}1.29& 6.12\\ 3.57& 8.46\\ 5.86& 7.60\\ 6.57& 8.30\\ 9.00& 6.90\end{array}\right],L_{5\&times;2}^{\left(old\right)}=\left[\begin{array}{cc}1& 1\\ 2& 5\\ 3& 3\\ 4& 4\\ 5& 2\end{array}\right],R_{5\&times;2}^{\left(old\right)}=\left[\begin{array}{cc}0.645& 0.500\\ 0.785& 0.070\\ 0.930& 0.710\\ 0.285& 0.860\\ 0.500& 0.430\end{array}\right],U_{5\&times;2}^{\left(old\right)}=\left[\begin{array}{cc}0.129& 0.100\\ 0.357& 0.814\\ 0.586& 0.542\\ 0.657& 0.772\\ 0.900& 0.286\end{array}\right]$

首先，用传统方法获得采样样本大小为7的扩展拉丁超立方采样。如图2，已有采样在新采样结构中的分布可知，竖列中有三列没有采样点，序号为2、4、6，横列中有三列没有采样点，序号为2、5、7，因此生成一个样本大小为3的LHS，采样矩阵为

$X_{3\&times;2}^{\left(1\right)}=\left[\begin{array}{cc}2.279& 2.547\\ 0.958& 1.965\\ 1.158& 0.971\end{array}\right]$

然后，将新生成的3个采样点映射到图2中，得到新采样点在新采样结构中的位置，累积概率值矩阵为

$U_{3\&times;2}^{\left(1\right)}=\left[\begin{array}{cc}0.817& 0.922\\ 0.149& 0.576\\ 0.549& 0.147\end{array}\right]$

进而，根据分布函数的反函数得到新采样点数值，如图5所示，新的采样矩阵为

$X_{(8\&times;2)}^{\left(new\right)}=\left[\begin{array}{cc}1.29& 6.12\\ 3.57& 8.46\\ 5.86& 7.60\\ 6.57& 8.30\\ 9.00& 6.90\\ 8.17& 9.01\\ 1.49& 7.70\\ 5.49& 6.36\end{array}\right]$

从图5可知，传统LHS扩展方法产生的新采样有8个点，不是严格意义上的采样大小为7的LHS。

下面采用本发明方法获得样本大小为7的扩展拉丁超立方采样，具体示意图如图1所示。

步骤1：计算已有采样点在新采样点结构中的分布。

计算映射到新采样结构中的排序矩阵

$L_{5\&times;2}^{\left(1\right)}=\left[\begin{array}{cc}1& 1\\ 3& 6\\ 5& 4\\ 5& 6\\ 7& 3\end{array}\right]$

b.计算映射到新采样结构中的随机数矩阵

$R_{(5\&times;2)}^{\left(1\right)}=\left[\begin{array}{cc}0.903& 0.700\\ 0.499& 0.698\\ 0.102& 0.794\\ 0.599& 0.404\\ 0.300& 0.002\end{array}\right]$

步骤2：计算需要保留的采样点矩阵

根据构建采样点间的邻接矩阵A_5×5，将采样点间的关系用矩阵中的元素表示；

$A_{5\&times;5}=\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 1& 0\\ 0& 1& 1& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$

从1到5的所有整数构成排列组合矩阵记为B_p，B_p的每一行是一个排列组合,则

$B_p={\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0\\ 2& 0& 0& 0& 0\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ 1& 2& 3& 4& 5\end{array}\right]}_{120\&times;5};$

根据B_P中排列组合的序号剔除对应采样点，即邻接矩阵A_5×5中对应的行和列，得到新的邻接矩阵A'；当第五次时，剔除第四个点之后得到邻接矩阵为零矩阵。则非零元素构成的向量B'_p＝[1235]，B'_p中的元素就是剔除采样的序号。A'的维数n'＝4，即保留n'个采样点；

c.按照向量B'_p元素的值剔除和中对应的行，分别得到

$L_{4\&times;2}^{\left(2\right)}=\left[\begin{array}{cc}1& 1\\ 3& 6\\ 5& 4\\ 7& 3\end{array}\right]R_{4\&times;2}^{\left(2\right)}=\left[\begin{array}{cc}0.903& 0.700\\ 0.499& 0.698\\ 0.102& 0.794\\ 0.300& 0.002\end{array}\right]U_{4\&times;2}^{\left(2\right)}=\left[\begin{array}{cc}0.129& 0.100\\ 0.357& 0.814\\ 0.586& 0.542\\ 0.900& 0.286\end{array}\right]X_{4\&times;2}^{\left(2\right)}=\left[\begin{array}{cc}1.29& 6.12\\ 3.57& 8.46\\ 5.86& 7.60\\ 9.00& 6.90\end{array}\right];$

步骤3：计算新增采样点矩阵即

生成(n+m)×r维的新排序矩阵L'_7×2，其中L'_7×2的每一列都是随机生成的整数1到7的一个排列, $L_{7\&times;2}^{\&prime;}={\left[\begin{array}{ccccccc}6& 3& 7& 5& 1& 2& 4\\ 4& 7& 1& 6& 2& 3& 5\end{array}\right]}^T;$

将L'_7×2每一列中与相应列相同的元素删除，剩余元素构成的矩阵记为

$L_{3\&times;2}^{\left(3\right)}=\left[\begin{array}{cc}6& 7\\ 2& 2\\ 4& 5\end{array}\right];$

生成随机数矩阵

$R_{3\&times;2}^{\left(3\right)}=\left[\begin{array}{cc}0.800& 0.916\\ 0.649& 0.734\\ 0.743& 0.392\end{array}\right];$

计算累积概率值矩阵

$U_{3\&times;2}^{\left(3\right)}=\frac{L_{3\&times;2}^{\left(3\right)}-\stackrel{\&RightArrow;}{1}+R_{3\&times;2}^{\left(3\right)}}{7}=\left[\begin{array}{cc}0.829& 0.988\\ 0.240& 0.248\\ 0.535& 0.628\end{array}\right];$

由公式求得采样值矩阵

$X_{3\&times;2}^{\left(3\right)}=\left[\begin{array}{cc}8.29& 9.61\\ 2.40& 6.76\\ 5.35& 7.84\end{array}\right].$

步骤4：计算ELHS采样矩阵

$X_{7\&times;2}^{\left(new\right)}=\left[\begin{array}{c}{X}_{4\&times;2}^{\left(2\right)}\\ X_{3\&times;2}^{\left(3\right)}\end{array}\right]={\left[\begin{array}{ccccccc}1.29& 3.57& 5.86& 9.00& 8.29& 2.40& 5.35\\ 6.12& 8.46& 7.60& 6.90& 9.61& 6.76& 7.84\end{array}\right]}^T,$

得到采样点为7的扩展拉丁超立方采样，采样结果见附图6。

本发明提出的扩展LHS方法可以应用于多种领域，例如电力系统概率潮流计算、混凝土抗压强度预测模型建立和航空器性能的灵敏度分析等。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。

Claims (4)

1.一种拉丁超立方采样的扩展方法，其特征在于，包括以下步骤：

步骤1：计算已有采样点在新采样结构中的分布，包括计算原始排序矩阵映射到新采样结构中的排序矩阵原始随机数矩阵映射到新采样结构中的随机数矩阵其中排序矩阵的矩阵元素其中i＝1,2,…,n，j＝1,2,…,r，为向下取整；随机数矩阵的矩阵元素其中i＝1,2,…,n，j＝1,2,…,r，n为已有样本大小，n+m为扩展得到的样本大小，n、m、r均为自然数；

步骤2：计算需要保留的采样点矩阵首先根据构建采样点间的邻接矩阵A_n×n，将采样点间的关系用矩阵中的元素表示；然后对采样点1,2,…,n的每个组合进行全排列，全排列矩阵为B＝[B₁B₂…B_q]^T；取p＝1，根据B_P中排列组合的序号剔除对应采样点，即A_n×n中对应的行和列，得到新的邻接矩阵A'，如果A'不为零阵，则p＝p+1，重新删除采样点构建；如果A'为零矩阵，记B'_p为B_P中非零元素构成的向量，则B'_p中的元素就是剔除采样的序号，A'的维数记为n'，即保留n'个采样点；最后按照向量B'_p元素的值剔除和中对应的行，分别得到和其中邻接矩阵A_n×n的矩阵元素为

步骤3：计算新增采样点矩阵首先生成(n+m)×r维的新排序矩阵L'_(n+m)×r，将L'_(n+m)×r的每一列与相应列中相同的元素删除，将剩余元素构成的矩阵记为然后生成随机数矩阵由此计算累积概率值矩阵最后求得新增采样值矩阵

步骤4：计算ELHS采样矩阵根据步骤1、2、3得到新的采样值矩阵即为样本大小为(n+m)的LHS，其中采样值矩阵 $X_{(n+m)\&times;r}^{\left(new\right)}=\left[\begin{array}{c}{X}_{{n_{\&prime;}}^{\&prime;}\&times;r}^{\left(2\right)}\\ X_{(n+m-n^{\&prime;})\&times;r}^{\left(3\right)}\end{array}\right],$ 其中分别由步骤2、步骤3求得。

2.根据权利要求1所述的拉丁超立方采样的扩展方法，其特征在于，步骤2中所述的全排列矩阵为B＝[B₁B₂…B_q]^T，其中B的每一行都是1,2,…,n的一个排列，不足n个元素的最后用0补齐，且B_i中0元素的个数不少于B_j，1≤i≤j≤n。

3.根据权利要求1或2所述的拉丁超立方采样的扩展方法，其特征在于，步骤3中所述的新排序矩阵L'_(n+m)×r的每一列都是随机生成的整数1到n+m的一个排列，随机数矩阵的矩阵元素r_i,j为[01]上均匀分布的随机数，累积概率值矩阵的矩阵元素的计算公式为

$u_{i,j}^{\left(2\right)}=\frac{l_{i,j}^{\left(2\right)}-1+r_{i,j}^{\left(2\right)}}{(n+m)}.$

4.根据权利要求1或3所述的拉丁超立方采样的扩展方法，其特征在于，步骤3中所述的新增采样点矩阵的计算公式为x_i＝F_i ^-1(Prob_i)。