使用阻尼BFGS更新构建Metropolis Hastings提案

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

英文标题：
《Constructing Metropolis-Hastings proposals using damped BFGS updates》
---
作者：
Johan Dahlin, Adrian Wills, Brett Ninness
---
最新提交年份：
2018
---
英文摘要：
  The computation of Bayesian estimates of system parameters and functions of them on the basis of observed system performance data is a common problem within system identification. This is a previously studied issue where stochastic simulation approaches have been examined using the popular Metropolis--Hastings (MH) algorithm. This prior study has identified a recognised difficulty of tuning the {proposal distribution so that the MH method provides realisations with sufficient mixing to deliver efficient convergence. This paper proposes and empirically examines a method of tuning the proposal using ideas borrowed from the numerical optimisation literature around efficient computation of Hessians so that gradient and curvature information of the target posterior can be incorporated in the proposal.
---
中文摘要：
基于观测到的系统性能数据计算系统参数及其功能的贝叶斯估计是系统辨识中的一个常见问题。这是之前研究的一个问题，其中使用流行的Metropolis-Hastings（MH）算法对随机模拟方法进行了研究。之前的这项研究发现了调整{建议分布，因此MH方法提供了充分混合的实现，以提供有效的收敛。本文提出并实证检验了一种方法，使用从数值优化文献中借用的有关Hessians有效计算的思想来调整建议，以便将目标后验的梯度和曲率信息纳入提议
---
分类信息：

一级分类：Statistics        统计学
二级分类：Computation        计算
分类描述：Algorithms, Simulation, Visualization
算法、模拟、可视化
--
一级分类：Quantitative Finance        数量金融学
二级分类：Computational Finance        计算金融学
分类描述：Computational methods, including Monte Carlo, PDE, lattice and other numerical methods with applications to financial modeling
计算方法，包括蒙特卡罗，偏微分方程，格子和其他数值方法，并应用于金融建模
--

---
PDF下载：
--> Constructing_Metropolis-Hastings_proposals_using_damped_BFGS_updates.pdf (630.2 KB)

2022-6-2 20:46:19 上传

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Metropolis Polis ASTIN Metro Etro

使用阻尼BFGS构建大都会黑斯廷斯提案更新Johan Dahlin、Adrian Wills和Brett Ninness*2018年5月9日摘要基于观测到的系统性能数据计算系统参数及其功能的贝叶斯估计是系统识别中的常见问题。这是一个先前研究过的问题，其中使用流行的Metropolis-Hastings（MH）算法对随机模拟方法进行了研究。之前的这项研究已经确定了一个公认的调整提案分布的困难，因此MH方法提供了高效混合的实现，以实现高效的收敛。本文提出并实证检验了一种调整方案的方法，该方法使用了从数值优化文献中借用的关于黑森人高效计算的思想，以便将目标后方的梯度和曲率信息纳入方案中。关键词：贝叶斯参数推理、状态空间模型、拟牛顿、BFGS。*给作者的电子邮件地址：firstname。lastname@newcastle.edu.au.JD和AW就读于澳大利亚纽卡斯尔大学工程学院。BN在澳大利亚纽卡斯特大学工程与建筑环境学院工作。这项工作得到了澳大利亚研究委员会发现项目DP140104350.1简介的支持。状态空间模型（SSM）在许多科学学科中普遍存在，包括系统识别【Ljung，1999】和金融【Durbin和Koopman，2012】。SSMs中的一个常见问题是参数θ的估计∈ Rnθ给出了一些观测数据y1:T={y，···，yT}。在本文中，我们考虑用x表示的ssm的这个问题~ uθ（x），xt+1 | xt~ fθ（xt+1 | xt），yt | xt~ gθ（yt | xt），（1），其中|θ、fθ和gθ表示由θ参数化的已知密度。这里，xt∈ RNX和yt∈ rNude记录时间t时系统的状态和观测值。

请注意，此参数化（1）包括大多数非线性和非高斯SSM以及输入ut∈ rnuca可以作为fθ和gθ的参数添加。估算（1）中θ的一种常用方法是采用隐含的一步预测分布pθ（yt | yt-1） 形成似然π（θ）=pθ（y1:T）=pθ（y）TYt=2pθ（yt | y1:T-1） ，（2）的观测数据。然后，通过似然的最大化参数给出θ的估计，即bθML=arg maxθpθ（y1:T），（3），这是众所周知的最大似然（ML）估计【Ljung，1999】。实际上，pθ（yt | yt-1） 可用于大多数SSM，但可使用粒子方法进行无偏估计【Doucet和Johansen，2011年】，这些方法可用于基于梯度的优化【Poyiadjis等人，2011年】和期望-最大化算法【Sch¨on等人，2011年】中，以近似解决（3）中的非对流可处理问题。在本文中，我们采用另一种方法，通过使用贝叶斯范式来估计θ【Peterka，1981，Robert，2007】。这相当于计算后验概率（θ| y1:T）∝ pθ（y1:T）p（θ），（4），其中p（θ）是θ的先验分布，可用于编码有关参数的先验用户信息。然而，这种后验概率很难计算，因为无法以闭合形式计算可能性（2）。相反，我们利用一种随机模拟方法来解决这个难题，它构造了一个随机数生成器，使得θk~ π(θ). （5） 这是统计文献中广泛应用的方法，在过去十年中，它在应用中的使用呈爆炸式增长。构建合适的随机数生成器的标准方法是采用（粒子）Metropolis-Hastings（MH；Robert and Casella，2004，Andrieu et al.，2010），这是一种非常通用的算法，用于计算π（θ）的实现，前提是可以逐点计算。

不幸的是，实现合理的收敛需要仔细调整算法，这本质上需要利用未知后验值的信息。这一困难在文献中得到了很好的认可，通常通过采用自适应方法【Andrieu和Thoms，2008】和几何信息的加入【Girolami和Calderhead，2011】来缓解。后一种方法需要计算对数后验的Hessian值，即使对于适用标准Kalman方法的线性高斯SSM，也很难直接计算。当使用粒子方法【Doucet和Johansen，2011】来估计潜在状态和可能性时，问题就更糟了。这是经验观察的结果，即使用粒子方法获得的Hessian估计值通常是有噪声的，甚至使用大量粒子也不准确。然而，有时即使使用少量粒子，梯度估计也很精确。因此，研究利用噪声梯度信息估计黑森值的问题很有意义。本文的贡献是探索使用阻尼Brodyen–Fletcher–Goldfarb–Shanno（BFGS）更新，仅使用梯度信息近似由Hessian编码的局部曲率。Ths是解决光滑数值优化问题的一种广泛使用的方法【Nocedal和Wright，2006年】。有关在MH内使用BFGS的相关工作包括Zhang和Sutton【2011年】，作者将这一想法应用于回归问题。Dahlin等人【2015b】将这项工作扩展到一类具有难处理可能性的SSM。本文的主要创新之处在于使用阻尼BFGS更新，以确保即使采用粒子方法，Hessian也是正半有限的。

最后，与早期在MH内使用BFG的尝试相比，所提出的方法提供了优异的性能。因此，我们获得了良好的MH性能，无需繁琐的用户调整，这是迈向自动贝叶斯推理方法的一步。2后验抽样存在一套所谓的马尔可夫链蒙特卡罗（MCMC；Robert和Casella，2004）方法，用于构建马尔可夫链，生成具有用户指定不变分布π（θ）的实现{θk}。由于在温和的假设下，马尔可夫链的实现具有收敛于链的变分分布的分布，因此这提供了一种构建具有任意目标分布π的随机数生成器（5）的方法。然后可以使用这些样本来近似后验值。给出后验点估计，如条件平均值bθCM=Eθy1:T=可以得到Zθπ（θ）dθ，（6）。这些都是有趣的，因为它们具有最小均方误差特性，并且不依赖于作为ML估计的渐近结果（3）。此外，通过计算（边缘）后验密度（θi | y1:T）=Zπ（θ）dθ，可以获得估计参数向量的每个元素i的误差界-i、 （7）式中θ-IDE注意到没有第i个元素的向量θ。不幸的是，（6）和（7）都需要对多维积分进行评估，这在计算上可能具有挑战性，尤其是当nθ=dim{θ}增长时。这导致任何任意（可测量）函数的期望值Д：Rnθ→ 由π[Д]=E[Д（θ）]=ZД（θ）π（θ）dθ给出的R，可以使用随机数生成器中的样本用bπK[Д]=KKXk=1Д（θK），（8）来近似。此外，我们还发现，估计量服从强大数定律，并且是一致的，即bπK[Д]a.s。-→ π[Д]，K→ ∞.

（9） 选择Д（θ）=θ，然后给出条件平均值估计（6）的近似值，选择Д（θ）作为适当的指示函数，然后给出π（θ）的近似值作为样本直方图。MH算法可以说是实现算法1 Metropolis Hastings（MH）输入：K>0，θ和q的最广泛使用的MCMC技术之一。输出：{θ，…，θK}。1： 计算π（θ）。2： 对于k=1至k do3：样品θ~ q（θ|θk-1） 使用（10）。4： 使用卡尔曼或粒子方法计算π（θ）。5： [0，1]上的ωkuniformly样本。6： 如果ωk≤ min{1，α（θ，θk-1） }由（11）给出。然后7：接受θ，即θk← θ.8： else9：拒绝θ，即θk← θk-1.10：结束if11：结束（5）。它通过采用任意建议马尔可夫链q（θk |θk）进行操作-1） 并通过随机接受来自该基链的实现来调节它。接受概率取决于目标π（θ）。在理论上（但不是在实践中），可以非常自由地选择提案分布，但在大多数情况下，使用高斯提案，qθ|θk-1.= Nθ; uθk-1., Σθk-1., （10） 其中θ表示候选参数，u和∑分别表示均值和协方差函数。通过（10）生成候选参数后，接受概率为α（θ，θk-1) = 1 ∧π（θ）π（θk-1） q（θk-1 |θ）q（θ|θk-1） ，（11）其中∧ b=最小值（a，b）。我们设置θk← θ如果候选参数被接受，θk← θk-1如果弹出。注意，我们只需要能够逐点计算π（θ）就可以实现算法1中所述的MH。一个关键点是，收敛速度（9）取决于实现{Д（θk）}的相关性。不相关程度越高，收敛速度越快，因此对于给定的有限个实现数K，近似值（8）越好。

这在MCMC文献中得到了很好的理解，其中随机近似（8）中的方差与积分自相关f=1+2成正比∞Xk=2corr{Д（θ），Д（θk）}，（12），也称为效率因子（IF）。反过来，实现的相关性{Д（θk）}关键取决于方案q（θ|θk）的选择-1).从Girolami和Calderhead【2011】的工作中可以看出，通过在提案（10）中加入关于后验概率的梯度和曲率信息，可以极大地改善混合（即自相关）。当后验面非各向同性时，这一点尤为重要，即某些参数对后验面值的影响程度大于其他参数。这种影响也会影响参数空间，这使得曲率的局部信息对于增加混合非常重要。3 Hessian估计一个重要的进一步问题是，如引言所述，当MH内采用粒子方法时，很难以有效的方式获得一般SSM的曲率信息。对数后验梯度很容易使用Fisher恒等式进行有效估计【Capp\'e等人，2005年】，G（θ）= 对数pθ（y1:T）= Eθhlog pθ（x1:T，y1:T）y1:Ti，（13），其中可以使用卡尔曼或粒子平滑器来计算或近似期望值。Louis’identity[Capp'e等人，2005年]可以以类似的方式用于估计对数后验的负hessian。不幸的是，由此产生的估计器通常会受到较大噪声敏感性的影响，这会导致正不确定性的频繁损失。这是有问题的，因为hessian的倒数经常作为协方差包含在提案（10）中。另一种方法是通过试运行计算Hessian估计值。

这将产生一个所谓的预处理矩阵，可用于缩放提案。然而，Girolami和Calderhead【2011年】、Nemeth等人【2016年】和Dahlin等人【2015a】等的数值研究结果表明，这种与参数无关的方法在混合方面是次优的。3.1阻尼有限内存BFG为了解决这些问题，我们建议利用优化文献中的知识，其中曲率信息通常是根据梯度信息估计的。这是非常成功的拟牛顿算法中使用的方法，该算法允许优化非线性函数。同样，与后验的Hessian估计相比，使用粒子方法获得的梯度估计通常是准确的，并且计算成本相对较低。在经典数学优化中，开发了一类所谓的拟牛顿方法【Nocedal and Wright，2006】，将曲率信息纳入搜索方向计算中，以加速收敛。这些方法的细节可以在标准参考文献中很容易找到，如Nocedal和Wright【2006】，但本质上，这些方法利用梯度和迭代信息来形成Hessian或其逆的估计。最著名的拟牛顿方法之一是BFGS方法，该方法对当前的Hessian估计Hl进行秩2更新，通过递归Hl+1=（I）形成更好的估计Hl+1- ρlslz>l）Hl（I- ρlzls>l）+ρlzlz>l，（14）ρl=（z>lsl）-1，sl=θl- θl-1，zl=G（θl）- G（θl-1).可以观察到，如果所有迭代的ρl>0，Hessian估计将保持正定义。在优化中，通过满足Wolfe条件的a线搜索算法可以保证将该条件设置为ρlca，参见【Nocedal和Wright，2006年，第8章】。

不幸的是，在MH设置中，这无法实施，因为这将导致马尔可夫链收敛到单个点。为了改善这个问题，我们在这里采用了所谓的阻尼BFGS方法，其中Zli被rlviarl=βlzl+（1- βl）Hlsl，βl=1，如果s>lzl≥ 0.2s>lHlsl（0.8s>lHlsl）/（s>lHlsl- s> lzl），如果s>lzl<0.2s>lHlsl，Hl+1=（I- ρlslr>l）Hl（I- ρlrls>l）+ρlrlr>l。除阻尼项外，我们还采用有限内存【Nocedal和Wright，2006，第9章】实现阻尼BFGS方法，以便计算负载保持适度。4拟牛顿提案要构建一个好的MH提案，将在输入（10）的均值和协方差函数中使用梯度信息以及BFGS的Hessian估计。对数后验曲线的二阶泰勒展开（Dahlin et al.，2015a）得出的一个典型选择是u（θ）=θ+H-1（θ）G（θ），∑（θ）=H-1（θ），（15）式中 > 0表示用户指定的步长。因此，我们可以将该建议视为后验值的局部高斯近似，这应该允许从中进行有效采样。激励（15）的另一种方式是将其视为黎曼流形上的随机游动，参见Girolami和CalderheadAlgorithm 2拟牛顿命题：ψk，M，{θi，G（θi）}ki=k-Mandδ>0。输出：θ。1： 从ψk中提取\'M唯一元素，并按升序（相对于对数目标）对其排序，以获得\'ψk，M.2：如果\'M≥ 2然后3：初始化Hessian估计值H.4：对于l=1到'M do5：根据'ψk，M中的lth对计算sland zl。6：执行更新（14）以获得Hl。7： 结束8：Set∑QN（ψk，M）=-H'M（θ）。9： else10：设置∑QN（ψk，M）=δIp。11： 结束if12：从（16）中取样以获得θ。详情参见【2011年】。将BFGS算法用于估计Hessian需要我们对toMH进行一些更改。

主要问题是，算法中来自M次迭代的信息被用于构建提案分布。在MH的标准版本中，由于马尔可夫特性，建议中只允许使用上一次迭代的信息。为了解决这个问题，我们需要将马尔可夫链从一阶链扩展到anM阶链。这允许MH保留其有效性，如Zhang和Sutton【2011】以及Dahlin等人【2015b】所述。算法1的主要算法变化是在步骤4中使用平滑器和算法2计算梯度和Hessian。此外，MH中的建议步骤由Q中的抽样代替θ|ψk，M= Nθ; uQNθk-M, Σ-1QNψk，M, （16） uQN（θk-M） =θk-M级+Σ-1QN（ψk，M）G（θk-M） 。使用算法2中的过程，ψk，M，{θi，G（θi）}ki=k-M、 最后，我们将算法1中的步骤9更改为θk← θk-M当候选参数因提案现在集中在θk上而被拒绝时-M、 5数值说明我们提供了三个数值说明，以了解所提出的算法，并将其与文献中的其他备选方案进行比较。附录A总结了实现细节，可以按照第6节所述下载源代码。时间标签。规则。如果Iter符合Cor.max。