SMC-ABC方法估计随机模拟模型

在案件中→ 0，所以Kh（ty-t） 在原点（即ty=t）成为点质量，在其他地方为零，如果ty对θ有效，那么难以处理的后边缘πM（θ| ty）=π（θ| ty）=π（θ| y）将被精确地恢复（尽管小h通常是不切实际的）。文献中研究了各种平滑内核的选择（例如Marjoram等人[2003]；Beaumont等人[2002]；Peters等人[2012]）。在我们讨论无似然抽样法或ABC抽样法时，可以方便地考虑联合分布（16）的推广≥ 1辅助汇总向量πJ（θ，t1:S | ty）∝~Kh（ty，t1:S）f（t1:S |θ）π（θ），其中t1:S=（t，…，tS）和t，tS~ f（t |θ）是由（难处理的）模型生成的独立数据集。由于辅助数据集通过构造是条件独立的，给定θ，我们有f（t1:S |θ）=QSs=1f（ts |θ）。我们遵循Del Moral等人[2012]的方法，将内核K指定为~Kh（ty，t1:S）=S-1PSs=1Kh（ty- 产生关节后部πJ（θ，t1:S | ty）=cJ“SSXs=1Kh（ty- ts）#“SYs=1f（ts |θ）#π（θ），（18）cJ>0适当的归一化常数，其中在（18）中，我们扩展了K（ty）的一致核选择- ts）由Del Moral等人[2012]在一般情况下提出。很容易看出，通过构造，RTSπJ（θ，t1:S | ty）dt1:S=πM（θ| ty）将分布（17）视为边际分布（c.f.Del Moral等人[2012]）。S=1，πJ（θ，t1:S | ty）=π（θ，t | ty）对应于无似然设置中更常见的后关节（16）。对于πM（θ| ty）的后验模拟，有两种明显的方法≈ π（θ| ty）作为π（θ| y）的近似值。第一种方法是直接在增强模型πJ（θ，t1:S | ty）上采样，实现联合采样（θ，t1:S）∈ 在t1:S上进行后验边缘化之前（即，通过从采样器输出中丢弃TSA实现）。

在这种方法中，求和量t1:sar被视为增广模型中的参数。第二种方法是直接从低维空间πM（θ| ty）采样，通过蒙特卡罗积分近似积分（17），代替πM（θ| ty）的每个后验评估。在这种情况下，πM（θ| ty）∝ π（θ）ZTKh（ty- t） f（t |θ）dt≈π（θ）SSXs=1Kh（ty- ts）：=πM（θ| ty），（19）式中t，tS~ f（t |θ）。该表达式由多位作者（如Marjoram等人[2003]；Reeves and Pettitt[2005]；Ratmann等人[2009]；Toni等人[2009]；Peters等人[2012]）研究，需要多个生成的数据集t，tS，对于边际后验分布的每次评估πM（θ| ty）。与标准蒙特卡罗近似法一样，Var[^πM（θ| ty）]随着LIM的增加而减小→∞Var[^πM（θ| ty）]=0。对于边际后验分布，数量t1:s仅作为估计πM（θ| ty）的一种手段，否则不会明确进入模型。样本数量直接影响估计的方差。3通过基于种群的采样器估计贝叶斯LOB随机代理模型引入基于种群的无似然采样器，以避免MCMC采样器中的不良混合（Sisson等人[2007]；Toni等人[2009]；Beaumont等人[2009]；Peters等人[2012]；Del Moral等人[2012]）。这些采样器传播粒子群θ（1），θ（N），通过一系列相关密度φ（θ），φT（θT），定义了从分布φ到目标分布φT的平滑过渡，从分布φ可以直接采样。对于无可能性或ABC取样器，φkis通过允许Khn（ty-t） 将更高的密度放置在≈ t随着k的增加（也就是说，带宽hn随着n的增加而减少）。

因此，我们表示πJ，n（θ，t1:S | ty）∝~Khn（ty，t1:S）f（t1:S |θ）π（θ）和πM，n（θ| ty）∝ π（θ）RTS~Khn（ty，t1:S）f（t1:S |θ）dt1:Sfor n=1，T，分别在关节和边缘后部模型下。3.1近年来，基于序贯蒙特卡罗的SSMC方法已从工程、概率和统计学领域涌现出来。这些方法的变体有时出现在粒子过滤或相互作用粒子系统的名称下（例如，Ristic等人[2004]、Doucet等人[2001]、Del Moral[2004]），Crisan和Doucet[2002]、Del Moral[2004]、肖邦[2004]和K¨unsch[2005]对其理论性质进行了广泛研究。标准的SMC算法涉及找到一组滤波递归的数值解，例如非线性/非高斯状态空间模型产生的滤波问题。在这个框架下，SMC算法从分布πn（通常是自然发生的）序列中采样，n=1，T每个分布都定义在支持En=E×E×·······················································≥1每个定义的固定支持E。注：这不是一个产品空间，而是一个固定空间E。Del Moral等人[2006]、Peters等人[2005]、Peters等人[2012]和Targino等人[2015]将SMC算法推广到目标分布π都定义在同一支持E上的情况。这种推广称为SMC采样器，使SMC算法适应更流行的设置，其中状态空间E保持静态，也就是我们前面讨论的关于MCMC算法的设置。简而言之，SMC取样器从一系列分布πn（n=1，…）生成加权样本（称为粒子）。

，T，其中πT可能特别有趣。我们将目标分布称为πTas，例如模型参数的后验分布。程序上，从任意初始分布π中获得的粒子，以及一组相应的初始权重，通过三个过程依次通过序列中的每个分布π传播，包括突变（或移动）、校正（或重要性权重）和选择（或重采样）。分布为π皮重的最终加权粒子被视为来自目标分布π的加权样本。因此，给定一系列分布{πn（dθ）}Tn=1，目的是在nw（i）n，Θ（i）noNi=1表示的每个时间n，开发一个n加权随机样本的大集合，其中W（i）n>0，pni=1W（i）n=1。这些重要权重和样本，用nw（i）n，Θ（i）noNi=1表示，是已知的粒子（因此通常被称为粒子过滤器或相互作用粒子系统等算法）。为了使这些方法合理，我们需要通过这些样本构造的经验分布渐近收敛（N→∞) 每次n的目标分布π。这意味着对于任何πn可积分函数，例如用φ（θ）：E表示→ Rone将具有以下收敛性：NXi=1W（i）nφθ（i）na、 美国。→ Eπnφ(Θ). （20） 在SMC采样器算法中，开发了SMC算法的一种特殊变体，即SMC算法的一种改进。考虑由πn（θ），n=1，…，给出的一般分布序列，T，带θ∈ E、 其中，最终分布π是兴趣分布。通过引入后向核Lk的GA序列，得到了一个新的分布：1πn（θ，…，θn）=πn（θn）n-1Yk=1Lk（θk+1，θk）（21）可定义为粒子的路径（θ，…，θn）∈ 通过序列π，πn.对向后核的唯一限制是正确的边际分布是πn（θ）。

. . , θn）dθ，dθn-1=πn（θn）可用。在这个框架内，人们可以在标准SMC算法下处理构造的分布序列eπn。总之，SMC采样器算法包括三个阶段：1。突变，粒子从θn移动-1到θnvia突变核Mn（θn-1，θn）；2.修正，其中粒子相对于πnvia（增量重要性权重）重新加权（等式22）；3.选择，根据颗粒多样性的一些度量，通常是有效样本大小，可以对加权颗粒进行重新采样，以减少重要权重的可变性。更详细地说，假设在时间n-1，分布eπn-1可以用eπNn经验近似-1使用N加权粒子。这些粒子首先通过突变核Mn（θn）传播到下一个分布eπ-1，θn），然后分配新的权重Wn=Wn-1wn（θ，…θn），其中Wn-1是粒子在时间n时的重量-1和Wn是由Wn（θ，…，θn）=eπn（θ，…，θn）eπn给出的增量重要性权重-1（θ，…，θn）-1） Mn（θn）-1，θn）=πn（θn）Ln-1（θn，θn-1） πn-1（θn）-1） Mn（θn）-1，θn）。（22）产生的粒子现在是来自eπn的加权样本。因此，根据等式（22），在SMC取样器框架下，可以直接使用边缘分布πn（θn），例如wn（θ，…，θn）=wn（θn）-1，θn）。而反向内核的选择-1基本上，它们的规格会强烈影响算法的性能，这将在以下小节中讨论。因此，SMC采样器算法的基本版本将按照算法3.1中给出的方式明确进行。备注3.1在我们使用增量重要性抽样权重校正的所有情况下，表达式中的参数只需在标准化之前已知。

也就是说，完全可以接受的是，只能计算目标分布序列{πn}到正规化常数。只要所有粒子都存在相同的标准化常数，这是正确的，因为重正化步骤将纠正对重要性权重知识的缺乏。在实践中，这对此类方法的应用至关重要。顺序蒙特卡罗采样器1。初始化粒子系统；（a） 设置n=1；（b） 对于i=1，N、 绘制初始粒子Θ（i）~ p（θ）；（c） 评估增量重要性权重SNWΘ（一）取方程（22）并归一化权重，得到nw（i）o。按{πt}Tn=2.2的顺序对每个分布重复以下步骤。如果有效采样大小（ESS）=PNi=1，则重新采样（a）w（i）n< NEFFI小于阈值Neff，然后通过多项式或分层方法通过加权样本的经验分布对粒子进行重新采样；参见K¨unsch[2005]和Del Moral[2004]关于无偏重采样方案的讨论。突变和校正（a）设置n=n+1，如果n=T+1，则停止；（b） 对于i=1，N从变异核Θ（i）N中提取样本~ 锰Θ（i）n-1.;（c） 评估增量重要性权重Θ（i）n取方程（22）并归一化权重，得到nw（i）noviaW（i）n=W（i）n-1w（i）n（Θn）-1，Θn）PNj=1W（i）n-1w（i）n（Θn）-1，Θn）。（23）在这个阶段，对于希望使用这种SMC方法的从业者来说，更好地了解这类估计方法在此类数值近似精度方面的特性是有益的。我们简要介绍了一些已知的结果，这些结果基于最近的粒子方法浓度不等式示例，这些示例是有限样本结果（见Del Moral等人的讨论和参考文献）。

[2013]).本文给出的指数集中不等式在一些一般状态空间中定义的粒子权重和变异核的正则条件下是满足的；见Del Moral[2004]关于这些情况的具体概率细节。利用平均场颗粒模型的浓度分析，可以获得以下指数估计值（见Del Moral[2004]的讨论）和其中的参考文献。注：在下文中，当使用分布或密度的N粒子近似值，如π时，我们将用πN表示它。定理3.1（有限样本指数浓度不等式）适用于任何x≥ 0，n≥ 0，任何人口规模N≥ 1、事件发生的概率isPrπNn（ψ）- πn（ν）≤cN1+x+√十、+C√N√十、≥ 1.-E-x、 （24）其中，N个粒子样本估计器的定义如下：πNn（~n）=NXi=1W（i）Nθ（i）n对于稳定的SMC算法，即对初始条件不敏感的算法，如我们前面讨论的，常数c和（c，c）不依赖于时间参数。我们还可以将目标分布的粒子估计和真实分布之间的差异限定为如下所示。考虑任何θ=（θi）1≤我≤还有吗(-∞, x] =Qdi=1(-∞, n=Rd中的θi]细胞，我们假设fn（x）=πn(-∞,x]FNn（x）=πNn(-∞,x].利用经验粒子构造的分布函数和分布序列{π，π。

，πT}，我们可以说明密度序列πT的分布函数的以下推论。推论3.1适用于任何推论≥ 0，n≥ 0，任何人口规模N≥ 1、以下事件发生的概率√NFNn- Fn≤ cpd（y+1）大于1- E-y、 这个浓度不等式确保了粒子重新分配函数FNTC收敛到Ft，几乎可以肯定是一致范数。3.2难处理似然贝叶斯模型的序贯蒙特卡罗采样器在ABC后验设置的SMC采样器的背景下对其进行规范化，粒子总体θn-1从分布φn中提取-1（θn）-1） 在时间n- 1被核Mn（θn）突变为φn（θn）-1，θn）。变异粒子θn的权重可以得到为Wn（θn）=Wn-1（θn）-1） wn（θn）-其中，对于边际模型序列πM，n（θn | ty），增量权重为wn（θn）-1，θn）=πM，n（θn | ty）Ln-1（θn，θn-1） πM，n-1（θn）-1 | ty）Mn（θn）-1，θn）≈πM，n（θn | ty）Ln-1（θn，θn-1） πM，n-1（θn）-1 | ty）Mn（θn）-1，θn），（26），其中，在（19）之后，将内核带宽设置为ABC容差级别我们得到πM，n（θn | ty）：=π（θ）SSXs=1Kn（ty）- ts）与基于S Monte Carlo drawst的πM，n（θn | ty）的（无偏）估计成正比，tS~ f（t |θn）。在这里-1（θn，θn-1） 是一个逆时核，描述粒子在n时刻从φn（θn）到φn的突变-1（θn）-1） 在时间n-1.与ABC-MCMC算法一样，增量权重（26）由“有偏”比^πM，n（θn | ty）/^πn组成-1（θM，n）-1 | ty）单位≥ 1.如果我们现在考虑联合模型πJ，n（θ，t1:S | ty）下的顺序蒙特卡罗采样器，自然变异核因子n[（θn-1，t1:Sn-1） ，（θn，t1:Sn）]=Mn（θn-1，θn）SYs=1f（tsn | ty）（与Ln类似-1） ，遵循（26）的形式，增量重量精确为Wn（θn）-1，t1:Sn-1） ，（θn，t1:Sn）=SPsKn（ty）- tsn）π（θn）Ln-1（θn，θn-1） SPsKN-1（ty）- tsn-1） π（θn）-1） Mn（θn）-1，θn）。

（27）因此，由于增量权重（26，27）是等效的，它们为边际模型和联合模型πM（θ| ty）和πJ（θ，t1:S | ty）引入了相同的SMC算法。因此，虽然以πM（θ| y）为目标的边缘采样器的应用在理论上有一定的偏差≥ 1，与之前一样，通过与关节空间目标πJ（θ，t1:S | ty）上的等效采样器关联，它们实际上是无偏的。我们注意到，理论上无偏取样器针对πM（θ| ty），适用于所有≥ 1，可通过仔细选择内核Ln获得-1（θn，θn-1). 例如，Peters[2005]、Peters等人[2012]和Targino等人[2015]都使用了byLn给出的次优近似最优核-1（θn，θn-1） =πM，n-1（θn）-1 | ty）Mn（θn）-1，θn）RπM，n-1（θn）-1 | ty）Mn（θn）-1，θn）dθn-1，（28），其中增量权重（26）近似为wn（θn-1，θn）=πM，n（θn | ty）/ZπM，n-1（θn）-1 | ty）Mn（θn）-1，θn）dθn-1.≈ πM，n（θn | ty）/NXi=1Wn-1（θ（i）n-1） Mn（θ（i）n-1，θn）。（29）在这种向后核选择下，现在所有S的权重计算都是无偏的≥ 1，因为近似^πM，n-不再需要（26）分母中的1（θ| y）。3.2.1自适应计划：通过退火容差计划选择ABC贝叶斯分布序列在本节中，我们考虑如何采用ABC后验分布，并构建SMC采样器所需的分布序列。也就是说，我们解决了如何开发ABC特定目标分布序列的问题。我们选择按照我们称之为ABC反向退火的方法来设计这个序列。特别地，我们构造了一个目标后验分布序列{φn}n≥0基于严格递减的公差值构造，通常由序列表示{n} n≥0以便> > . . . > n> …>T

我们通过考虑φn获得了ABC后验分布序列，φn与涉及核的ABC可能性有关。如果我们认为内核的带宽由K给出n（ty）- t） 然后，我们将逐步更加强调，即密度，即≈ t随着n的增加（即带宽n随n增加）。因此，在SMC取样器程序中，我们分别表示了在关节和边缘后部模型下可能考虑的两种ABC结构：πJ，n（θ，t1:S | ty）∝~Kn（ty，t1:S）f（t1:S |θ）π（θ）πM，n（θ| ty）∝ π（θ）ZTS~Kn（ty，t1:S）f（t1:S |θ）dt1:S（30）现在，这个过程使其自适应的方面是选择π和πn+1之间的差异大小，对于每个n∈ {1,2…，T}以及最终停止点。本文提出在每一步对ABC目标分布序列进行自适应。AIM是逐步在线选择一系列分布，以便下一个分布和上一个分布之间的差异由公差控制N序列由序列中先前目标分布的粒子近似的“精确度”或效率控制。一个好的近似值表示一个人可能会采取更大的步骤，而一个较差的近似值表示应该采取更小的步骤。在形式上，我们进行调整，以使新的公差n、 在迭代n中，通过基于粒子系统的分位数匹配程序生成。