响应面排序的序贯设计

我们比较了自适应Gap-SUR（3.2）和Gap-UCB方法（3.1）（γk=1）·√a）基于真实间隙的均匀间隙(·). 右图：通过（4.7）计算的K=400时的L（^C（K），C）在六种不同EI方法中的箱线图。表1还强调了在响应面之间进行区分的好处，因为Con18 Ruimeng Hu和Michael Ludkovski当前的M-Gap算法相对于Gap-SUR明显更差（损失高出约30%）。这些方法之间唯一的区别是，Gap SUR在200个样本中进行了146次采样，而同时进行的方法仅限于对每个响应进行100次采样。所有在整个X×L范围内进行优化的方法都将重点放在噪音更大的Y上，对其进行70–85%的采样（参见数据栏）。作为最后的比较，表1的最后一行报告了Gap-SUR方法在实际环境中的性能，其中还必须通过学习θi，s，σ来训练GP核K`。所有参数，包括被视为GP协方差结构金块的观测噪声σ，都是通过MLE估计的。由于训练会将额外的噪声引入固定的响应曲面，因此算法性能必然会降低，尤其是在交叉算法运行的变化方面。这可能表明静态GP模型在这里并不理想。表1还显示，经验EL（^C（K））和实际损失L（^C（K），C）指标是一致的，因此前者可以用作内部在线评估工具，以监控估计分类的准确性。这两个指标之间的不匹配是由模型规格错误造成的，因为不正确地推断出u（x）的协方差结构会导致过度乐观：EL<L。这个问题在很大程度上独立于抽样方案，更多地涉及建模框架，而不是toEI采集函数。4.4. 许多表面。

我们的下一个示例将处理一个更复杂的设置，其中L=5个曲面和一个二维输入空间X=[-响应参数（θ，θ，η，t`）u（x，x）2- 十、- 0.5x（4,6.5,23，-10） u（x，x）2（x- 1） +2倍- 2（7.5,7.5,475,60）u（x，x）2sin（2x）+2（1,8,2,1.9）u（x，x）8（x- 1） +8x- 3（8,8,8000,300）u（x，x）0.5（x+3）+16x- 6（8，4，2500，150）我们假设恒定的同调观测噪声`（x，x）~ N（0，σ`），σ`=0.5`. GPM模型具有可分离的各向异性5/2协方差函数和特定的超参数，以及固定的趋势t`。图4显示了相应的分类C。通过从每个样本（x，x）生成10个LHS样本，序列设计在K=50时初始化；在每个步骤中，使用随机化方法从LHS候选集T（大小=100）中选择采样位置-贪心法 = 0.1.图4左上角的面板显示了使用Gap SUR采集功能，在K=500个样本之后，估计的分类^C（K），其他面板显示了分配给每个`=1，5.如图所示，该算法在X×L上联合采样时具有很高的分辨力。在任何给定的分类边界上，该算法有效地只对五个响应中的两个进行了采样，内生地恢复了最佳与次最佳测试的概念。因此，来自Y\'的样品主要位于表面u\'和其他表面的边界附近。这些轮廓，在哪里`= u`- minj6=`uj=0，正是Gap EI指标所针对的区域。由于C和C的边界最长，因此选择了相对较多的样本（D=126，D=109）；相反，最小的一组是C，它只接收D=70个样本。表2显示了不同采集功能的相对性能。

具体而言，我们比较（i）均匀抽样；（ii）空间填充LHS采样；（iii）并行ΓBeststrategy（4.5）排序响应面的顺序设计19-2.-1 0 1 2-2.-1012x1x2（x）=1（x）=2（x）=3（x）=5（x）=4-2.-1 0 1 2-2.-1012x1x21x2123456（x）=1（x）=2（x）=3（x）=5（x）=4-2.-1 0 1 2-2.-1012x1x1x0246810（x）=1（x）=2（x）=3（x）=5（x）=4总体^C`=1`=2-2.-1 0 1 2-2.-1012x1x2012345（x）=1（x）=2（x）=3（x）=5（x）=4-2.-1 0 1 2-2.-1012x1x1x0246810（x）=1（x）=2（x）=3（x）=5（x）=4-2.-1 0 1 2-2.-图4。X上的二维排名=[-2, 2] × [-2，2]使用Gap SUR启发式。左上角的面板：实心黑线表示真实的C（x，x），红色虚线表示K=500时的估计等级C（K）（x，x）。其他面板显示了5个响应面的边缘设计（x，x）1:D`（K）。阴影表示估计的经验间隙B`（x，x），`=1，5.我们观察到，大多数样品被B`\' 0.实心曲线表示真实类别C（x，x）的边界。这类似于基于熵的抽样；（iv）Gap UCB和（v）Gap SUR。我们注意到，对于许多曲面，关键不一定是曲面之间的预算分配（这里，对于相同的σ`，最佳D`\'大致相等），而是有效地放置最适合每个曲面的样本位置。这种影响可以通过比较非适应性策略（即x和`）与并行的ΓBest策略（4.5）（针对分类边界，但在`）与Gap SUR/Gap UCB策略（针对不同指数的分类边界的不同部分）来观察。上述步骤中的每一步都会带来显著的性能提升；当观测噪声依赖于指数（或状态）时，它预计会更加明显。5.流行病管理案例研究。

最后一个例子是基于传染病背景下的控制问题[32,34,35,39]。考虑随机SIR模型，该模型是一个分区状态空间模型，将人口池划分为三类：易感人数St、感染人数St和恢复人数Rt。我们假设固定的人口规模EM=St+It+Rt，因此状态空间是二维单纯形X={（s，i）∈ Z+：s+i≤ M} 。在典型的环境中，M∈ [10,10]，因此X是离散的，但太大，无法明确列举20 Ruimeng Hu和Michael LudkovskiTable 2对于二维示例，Z（500）表示真实损失与经验损失。对于UCB启发法，冷却计划的形式为γk=c√log k.错误概率为ErrP rob=Ave（1-p（500）最佳（x））超过测试集。向量D`（500）列出了Y`，`=1，（4.64E-5）5.47E-3（2.39E-5）5.47E-3（2.39E-4）4.10（100100100100100100100）4.10（100100100100100100100）4.10（100100100100100100）4.10（100100100100100100100100100100100）4.10（100100100100100100100100100100）10）4.10（100100100100100100100100100100100100100100100100）10）10）10）4（100100100100100100100100100）10）10）4（100100）10）10）10）非自适应LHS 5（10）非自适应LHS 5.97E-3（10（10（10）非自适应LHS 5.97E-3（10（10（10）5（10（10）5（10）5.97E-3（10010010010010010010010010）非自适应LH5（10）5（10）5（10 Gap UCB，c=0.53.41E-3（1.45E-5）2.97E-3（1.14E-4）3.05%（129、103、104、72、92）（按|X |\'10的顺序）。（St，It）的动态是时间固定的，将在（5.3）中详细说明。控制员的目标是通过及时干预减轻疫情影响，例如通过降低个人接触率来降低感染率的社会距离措施；从数学上讲，这相当于修改（St，It）的动力学。

为了进行成本效益优化，我们一方面引入流行病成本，这里的成本与累积感染者的数量成正比，另一方面引入干预成本，这与当前CISt剩余易感者的数量成正比。然后，通过比较不采取行动的预期成本u（s，i）（以当前状态（s，i）为条件）和立即采取行动的预期成本uA（s，i），可以（短视地）优化干预方案。更准确地说，让u（s，i）：=E[s- ST |I=I，S=S]和（5.1）uA（S，I）：=EA[S- ST |I=I，S=S]+CIs。（5.2）以上，T=inf{T:It=0}是疫情的随机结束日期；由于固定的人群和感染后对疾病的假定免疫力，该流行病保证有一个婴儿期。不同之处在于- 因此，STTHES精确地测量了在暴发期间某个时间点感染的原始易感者的总数。总体目标是对u和uA进行排序，干预区域对应于{（s，i）：uA（s，i）>u（s，i）}。由于没有可用于u′s的分析公式，一个合理的程序（由于易于处理SIR模型的大量扩展，因此也是首选）是一个蒙特卡罗抽样器，给定初始条件s=s，I=I和区域∈ {0，A}生成一条轨迹（St，It）（ω），并使用它来计算路径St（ω），连接到（1.1）的框架。从政策角度来看，（5.1）-（5.2）中的权衡围绕着什么都不做，让疫情继续发展，这会为每个最终受感染的个体带来单位成本，或者实施预防性社会疏远措施，这会让每个易感者付出代价，但会降低未来受感染者的预期数量。典型的应对措施可能是公共广告宣传、学校停课或分发预防药物。

一般来说，一旦有足够大的疫情威胁，就需要进行干预。然而，如果干预成本较低，那么干预成本相对于其效益来说就太高了，因为疫情可能会自行结束。类似地，如果性病进展缓慢，易感人群自然会精疲力竭，再次使干预变得无关紧要（因为“为时已晚”）。量化这些场景需要一个精确的概率模型。（St，It）在各自规律下的动力学，并遵循连续时间马尔科夫序贯设计，对具有以下两个过渡通道的21条链进行排序：（感染：S+I）→ 2I，速率为βjStIt/M，j=0，A；康复：我→ R与它一起。）（5.3）如上所述，βA<β被解释为干预机制中受感染者和易感者之间的接触率降低，从而减少疫情增长和影响。（5.3）中描述的马尔可夫链（St，It）可以使用Gillespie时间步进算法[19]进行模拟，利用状态转换之间的逗留时间具有（依赖于状态的）指数分布，并且与下一个转换类型无关这一事实。然而，这些模拟非常耗时，需要O（M）均匀绘制。因此，有效的预期成本排序在应用中非常重要。备注5.1。由于（5.3）意味着每个感染期都有一个独立的Exp（γ）分布，因此E[S]- ST]=γERTItdt, 因此（5.1）也可以解释为与预期感染总天数的比例。我们注意到，在这个例子中，输入空间X是离散的，但这需要对算法1的实现进行最小的更改。最大的调整是，（1.2）中的噪声σ`（x）未知。σ`（x）’的知识对于训练GP协方差核K`至关重要，参见例（2.9）。

事实上，虽然可以同时训练K`和恒定观测噪声σ（后者在GP文献中被称为“金块”，可以通过最大似然法推断），但与状态相关的噪声K是无法识别的。我们通过分批程序（与[3，第3.1节]进行比较）来解决这个问题，以便随时估计σ`（x）。也就是说，我们使用相同的站点x≡ （s，i）r次，以获得独立样本y（1）`（x），y（r）`（x）来自相应的y`（x）。这允许估计条件方差eσ`（x）：=r-1rXi=1（y（i）（x）- y`（x）），其中y`（x）=rrXi=1y（i）`（x）是样本平均值。此外，如[40，第4.4.2节]所示，我们可以将x处的r样本视为具有噪声方差eσ`（x）/r的单一设计条目（x，\'y`（x））。由此产生的后平均设计尺寸减少了一个系数，大大加快了克里金模型拟合和更新的计算速度。形式上，算法1中的EI步骤被替换为使用（xk+1，`k+1）=（xk+2，`k+2）=（xk+r，`k+r）并每r次地面迭代重新计算一次EI分数。在我们的研究中，我们设定M=2000，β=0.75，βA=0.5，γ=0.5，每名易感者的干预成本为CI=0.25。图5显示了最终的决策边界C.在黑暗地区，干预的相对成本较低，因此行动是首选。例如，从atI=10，S=1800开始，如果不采取任何行动，疫情将影响40%以上的易感人群（预计成本约为800人），而在社交距离下，影响约为60名受感染者（导致60+CIS’510的总预期成本大大降低）。在光明地区，观望法的预期成本较低。例如，在I=50，S=1400时，未采取任何措施的新感染者的预期数量为385，而在0时，对抗措施的成本更高。25 × 1400 + 102 = 452.

总的来说，图5显示，最佳决策对当前易感者的数量非常敏感。这一特征是因为，当感染率主导恢复时（繁殖率R:=（β/γ）（S/M）高于1），疫情就会爆发。因此，对于一个敏感度超过85%（S>1700）的群体，初始增长率满足22 Ruimeng Hu和Michael Ludkovski1200 1300 1500 1600 1700 18000 50 100 150 200疑似感染（S，i）=0（S，i）=图5。拟合响应边界C使用Gap SUR预期改善指标的流行病响应示例。散点图表示K=200时的设计Z（K）；三角形表示初始设计Z（K），并围绕自适应放置的（s，i）K:K（绿色：Y；黄色：YA）旋转。βS/M>γ，可能引发暴发。然而，随着S的降低，接近βS/M\'γ的区域，随着疫情爆发的可能性和严重性的降低，社会距离变得不必要。特别是，图5显示，S<1350时，没有采取任何行动。在中间区域，有一个确定C（s，i）的非平凡分类边界。图5是通过使用Gap SUR采集功能构建自适应设计生成的，共有K=200个设计点，每个点的r=100批样品。输入空间被限制为X={s∈ {1200，…，1800}，我∈ {0, 200}}. 初始设计Z（K）包括50=25×2个位置，位于同一矩形5×5晶格上，每个Y，YA。在本例中，噪声水平σ`（s，i）高度依赖于状态，见图6。u表面的噪声要高得多，（s，i）\'（1800，5））的最大σ（s，i），而uA的最大噪声位于右上角。因此，Z（K）主要包含来自Yan的样本，并且在图的底部更密集。结论

在这篇文章中，我们构造了几个有效的顺序设计策略来解决L中最小值的确定问题≥ 2个响应面。我们的Gap SURheuristic将（1.1）与轮廓查找和贝叶斯优化相结合，为逐步减少不确定性框架提供了新的应用[10]。我们的Gap UCB启发式算法通过将X×L中所有可能的采样对视为ARM，并试图平衡armexploration和剥削来模拟多臂强盗。我们的方法基于克里格框架，但这主要是为了方便，并不重要。为此，可以使用最大似然法代替贝叶斯公式来拟合bu`（·），用点估计量及其标准误差代替后验M`（x）。因此，可以选择许多其他回归框架。然而，计算效率和顺序框架对建模u`（·）的可能方法造成了一些效率限制。一方面，我们需要很强的一致性，即各类别C（K）的收敛性→ C问→ ∞. 特别是，回归方法必须是非参数和局部的。另一方面，我们希望有一个顺序程序，允许在从^C（k）移动到^C（k+1）时高效地更新规则。最后，在实际环境中，进一步的挑战，如异方差、排序响应面的非高斯序列设计231200 1300 1400 1500 1700 1700 18000 50 100 150 200怀疑受影响1002003004005006001200 1300 1500 1700 1700 18000 50 100 150 200怀疑受影响1020304050Eσ（s，i）eσA（s，i）图6。在无对抗措施（左面板，`=0）和行动（右面板，`=A）制度下，流行病反应示例的估计噪声标准偏差eσ`（s，i）。

注意两个面板的不同色阶，带有σ（·） σA（·）适用于所有（s，i）。取样器Y〃和响应面的异质结构很重要。GP的一个合适替代方案是局部回归或黄土[44]，这是一个非参数回归框架，用于为u`（x）建立点态线性回归模型。黄土非常有效，非常适合具有未知噪声分布的异方差环境，如第5节所述。它还自动生成fit的后验均值和方差（允许使用基于bu`（x）和δ`（x）的衍生公式）。然而，黄土是不可更新的，如果要使用许多设计增强迭代，就会造成计算瓶颈。同时，安装速度非常快，因此根据实施情况，它可能仍然与更复杂的方法竞争。本着这种精神，分段线性回归（首先将X划分为几个单元，然后在每个单元中进行最小二乘回归）可通过Sherman Morrison Woodburyformulas进行更新，如果有明确的划分策略可用，则可采用分段线性回归。我们进一步注意到，GP克里金法只是构建实验设计的一种方便的临时替代方法。因此，一旦生成Z，人们可以切换到不同的响应面模型，以建立对u′s和^C的最终估计。例如，当基本平滑度（由协方差核指定）在X上强烈变化时，treed GP方法[24]允许响应面具有更高的精确性。因为treed GP模型的计算成本很高，在DoE期间使用vanilla GP和treed GP进行^C的最终估计，可以达成妥协。另一个富有成效的扩展是研究固定密度设置中的排序算法。