线性动态系统的噪声鲁棒在线推理

尽管存在许多其他重要的SMN分布（例如指数幂、对称α-稳定、逻辑、马蹄形、对称广义超对数分布）[34]，[35]，但在续集中，我们仅介绍上述特殊情况，因为相关的混合密度在计算上具有吸引力。在这种情况下，高斯分布可以被认为是一种退化的混合[33]DRAFT1）多元学生t分布：当X遵循一个多元学生t分布，其位置为u，标度为∑，自由度为ν（即X~ T（u，∑；ν）），X的pdf可以用以下SMN形式表示：p（X）=Z∞N（x |u，∑λ）Ga（λ|ν，ν）dλ，（7）其中Ga（·a，b）是形式Ga（λ| a，b）=baΓ（a）λa的伽马密度函数-1e-bλ，λ，a，b>0，（8）和Γ（a）是参数a>0的伽马函数。因此，X可以用以下分层形式表示X |（u，∑，ν，λ）~ N（u，λ）-1Σ), λ|ν ~ Ga（ν，ν）。（9） 2）Pearson VII型分布：如果X属于Pearson VII型家族，则相关密度由P（X |u，∑，ν，δ）=B（δ，）√νΣ1+（x- u)νΣ-（δ+1）/2，（10）其中u和∑是位置和比例参数，ν>0和δ>0是形状参数，b（a，b）是参数a>0和b>0的贝塔函数。皮尔逊VII型密度可以用asX |（u，∑，ν，δ，λ）分层表示~ N（u，λ）-1Σ), λ|(ν, δ) ~ Ga（ν，δ）。

（11） 当ν=δ时，我们回到学生的分布，如果ν=δ=1.3，则回到柯西分布：斜线分布可以按层次表示为x |（u，∑，ν）~ N（u，λ）-1Σ), λ|ν ~ Be（ν，1），（12）与0<λ<1和ν>0，其中Be（·）表示β分布。4） 多变量方差伽马分布：当X遵循位置为u、标度为∑、自由度为ν>0（即X）的多变量方差伽马（VG）分布时~ VG（u，∑；ν）），密度可以用以下分层形式x |（u，∑，ν，λ）表示~ N（u，λ）-1Σ), λ|ν ~ 其中IG（a，b）是由IG（λ| a，b）=baΓ（a）λ给出的反伽马密度-（a+1）e-b/λ。（14） 当ν=2时，VG变成拉普拉斯分布。德拉夫特。正态方差均值混合随机向量X遵循正态方差均值混合分布，如果可以表示为遵循[36]X=u+β∧+√∧Z，（15），其中∧独立于Z，Z根据协方差矩阵∑的零均值多元正态分布。随机变量∧是混合参数，分布在正实线和u，β上∈ R.以∧=λ为条件的X的分布是多元正态分布，给出为（X |∧=λ）~ N（x |u+βλ，∑λ）。请注意，当β=0时，NVMM将变成κ（λ）=λ的SMN。我们在下面描述了这类分布的一些特殊情况：1）广义双曲分布：如果∧按照广义逆高斯（GIG）分布asp（λ）=（a/b）p/22Kp分布，则（15）中的X遵循广义双曲（GH）分布(√ab）λp-1exp{-aλ+b/λ}，λ>0，（16），其中kp是第二类修正贝塞尔函数a，b∈ R+，和p∈ R.2）GH-skew-student的t分布：GH-skew-student的t分布的层次结构为[37]（X∧=λ）~ N（x |u+βλ，∑λ），λ|ν~ IG（ν/2，ν/2），（17），其中反伽马密度由（14）给出。

请注意，GH倾斜学生的t是GH分布的一个特例，其中GIG的参数选择为a=-ν/2（ν>0），b=√ν和p=0。此外，当β=0时，这将变成对称的Students t分布，如果ν→ ∞, 这变成了一个偏正态分布。3） GH方差伽马分布：GH方差伽马分布的层次结构为（X∧=λ）~ N（x |u+βλ，∑λ），λ|ν~ 其中伽马密度由（8）给出。备注1：注意，II-A和II-B中的混合参数∧是一个随机变量，∧的维数始终小于或等于随机向量X的维数。这一观察结果很重要，因为我们将在后面看到，PF可以用于在RBPF框架中瞄准低维∧。DRAFTIII。基于仿真的在线滤波器当动态系统（状态空间模型）是非线性和/或由非高斯噪声驱动时，滤波器密度p（xk | y1:k）通常在分析上难以处理。为了解决这个问题，随着时间的推移，人们提出了许多近似方法[38]。粒子滤波（PF）就是这样一种方法，它使用蒙特卡罗模拟来解决滤波问题。在本节中，我们将对PF和RBPF进行简要概述。A.粒子过滤（PF）在PF中，与密度p（x0:k | y1:k）相关的后验分布近似于一组N( 1） 加权粒子（样本）asbPN（dx0:k | y1:k）=NXi=1ew（i）kδx（i）0:k（dx0:k）；ew（i）k≥ 其中δx（i）0:k（A）是给定x（i）0:k和可测量集A的狄拉克度量，ew（i）kis是附加到每个粒子x（i）0:k的相关重量，因此pni=1ew（i）k=1。尽管分布bpn（dx0:k | y1:k）不允许关于勒贝格测度的定义良好的密度，但我们使用符号滥用来表示相关的经验密度asbpN（x0:k | y1:k）=NXi=1ew（i）kδ（x0:k）- x（i）0:k）。

（20） 虽然（20）在数学上并不严格，但从直觉上看，它比严格的测量理论符号更容易遵循，尤其是如果我们不关心理论收敛性研究的话。请注意，我们使用PF瞄准的后p（x0:k | y1:k）是未知的。（19）中的经验分布bpn（dx0:k | y1:k）是通过首先从建议分布π（x0:k | y1:k）生成样本x（i）0:kf获得的，然后使用重要性抽样的思想asw（i）k=p（x（i）0:k | y1:k）π（x（i）0:k | y1:k）获得相应的权重；ew（i）k=w（i）kPNj=1w（j）k.（21）给定此PF输出，现在可以近似地计算边际分布p（xk | y1:k）asbPN（dxk | y1:k）=PNi=1ew（i）kδx（i）k（dxk）。假设时间（k）-1） ，我们有后验概率（x0:k）的加权粒子近似-1 | y1:k-1） asbPN（dx0:k）-1 | y1:k-1） =PNi=1ew（i）k-1δx（i）0:k-1（dx0:k）-1). 现在有了一个新的量度k，我们希望用一组新的粒子（样本）来近似p（x0:k | y1:k）。标准PF使用以下后验路径空间递归P（x0:k | y1:k）∝ p（yk | x0:k，y1:k）-1） p（xk | x0:k）-1，y1:k-1） p（x0:k）-1 | y1:k-1). （22）对于马尔可夫状态空间模型，这变成了（x0:k | y1:k）∝ p（yk | xk）p（xk | xk）-1） p（x0:k）-1 | y1:k-1). （23）假设建议分布可以分解为π（x0:k | y1:k）=π（xk | x0:k）-1，y1:k）π（x0:k）-1 | y1:k-1） ，我们现在可以实现顺序重要性采样，通过从边缘提议核π（xk | x（i）0:k采样一个新的状态x（i）k，粒子被传播到时间kb-1，y1:k）和设置x（i）0:k，x（i）0:k-1，x（i）k. 随后使用（23），可通过w（i）k给出粒子的相应重量∝p（yk | x（i）k）p（x（i）k | x（i）k）-1） π（x（i）k |x（i）0:k-1，y1:k）ew（i）k-1.ew（i）k=w（i）kPNj=1w（j）k.（24）为了避免携带小重量的轨迹，并集中于大重量的轨迹，需要定期对粒子进行重新采样。

有效样本量Neff是一种衡量实际有助于近似分布的粒子数量的指标，通常用于决定何时重新采样。当NEFFR下降到指定阈值以下时，将执行重采样。文献中提出了许多有效的抽样方案。我们没有深入讨论细节，而是让感兴趣的读者参考[24]、[29]、[39]–[41]来了解PF的更全面的介绍。备注2：如果使用状态转移密度作为建议，并在每一步进行重采样，则相应的PF称为自举粒子滤波器。这很容易实现，因此在实践中非常普遍。B.Rao Blackwellized particle fi filtering（RBPF）虽然PF非常流行，而且已经存在了一段时间，但它在计算上要求很高，值得注意的是，它在扩展到更高维度时有严重的限制[30]。对于某些模型，当部分状态空间（有条件地）是可处理的时，就足以对状态空间的剩余不可控制部分使用PF。如果有可能利用这种分析子结构，则基于蒙特卡罗的估计将被限定在较低维度的空间。因此，得到的估计值通常比PF针对整个状态空间提供的估计值更好（就渐近方差而言），而且永远不会更差。由此产生的方法通常被称为Rao Blackwellized粒子过滤[32]，[42]–[45]。请注意，解析地（或使用解析近似法）求解状态向量的一部分，剩下的部分仅由PF作为目标。因此，RBPF可作为扩展到高维问题的实用工具。DRAFTIV。平稳非高斯环境下的LDS推断考虑（1）中的LDS，由已知的平稳非高斯噪声驱动，如第二节中所述，假设为分层高斯噪声。

给定该模型，并假设初始先验p（x）是已知的高斯密度，我们的目标是递归地针对难处理的滤波密度p（xk | y1:k）。设在时间步k，λwk和λek为wk和ek的相应混合参数。我们定义了辅助向量λk≡ （λwk，λek）T.然后噪声wk和ek是层次高斯aswk |λk~ Nuwk（λk），Qk（λk）, （25a）ek |λk~ Nuek（λk），Rk（λk）, （25b）其中uwk、uek是相应的平均值，qk和rkk是相应的分层高斯噪声协方差。这些参数可能取决于λk。这样的噪声表示允许CLGM，可以使用KF进行分析处理。这为在线推理问题提供了一个RBPF实现。在续集中，我们将描述这个RBPF框架。本节组织如下。我们首先概述了拟议RBPF框架的一个完整迭代，然后是相关PF针对的状态动力学（即混合参数）的规定。接下来，给出了该方法的算法描述。然后描述似然函数估计和p步超前预测。A.已知RBPf在时间零点p（x）处的后向传播周期。假设在时间步（k- 1） ，我们有联合目标分布p（xk-1，λ1:k-1 | y1:k-1). 这可以分解为aspxk-1，λ1:k-1.y1:k-1.= Pxk-1.λ1:k-1，y1:k-1.Pλ1:k-1 | y1:k-1., （26）式中p（λ1:k）-1 | y1:k-1） 以PF为目标，由一组N（>>1）加权随机粒子asp（λ1:k）经验给出-1 | y1:k-1） =NXi=1ew（i）k-1δ（λ1:k）-1.- λ（i）1:k-1） ，（27）更一般地，p（x）可以被认为是一个层次高斯密度。如果其中一个噪声是高斯噪声，则随后描述的推理框架仍然有效，而不是两个噪声都是非高斯噪声。

在这种情况下，λk表示一个辅助随机变量，代表非高斯噪声的混合参数。带ew（i）k的Draft-1.≥ 0和Pni=1ew（i）k-1= 1. 利用（27），我们得到了p（λk）-1 | y1:k-1） =NXi=1ew（i）k-1δ（λk）-1.- λ（i）k-1). （28）注意到，给定一个序列λ（i）1:k-1.动态系统现在成为CLGM。所以p（xk-1 |λ（i）1:k-1，y1:k-1） 可以通过KF获得，给出byp（xk-1 |λ（i）1:k-1，y1:k-1） =Nbxk-1（i），Pk-1（i）, （29）其中，为了符号清晰，我们抑制了参数对λ（i）1:k的依赖性-1.由于KF为每个序列运行（此后由索引i表示），因此在任何给定时间，共有N个KF并行运行。现在将（27）和（29）一起使用，滤波器分布p（xk-1 | y1:k-1） 可获得asp（xk-1 | y1:k-1） =Zpxk-1.λ1:k-1，y1:k-1.×pλ1:k-1.y1:k-1.dλ1:k-1.≈NXi=1ew（i）k-1Nbxk-1（i），Pk-1（i）, （30）是N个高斯分布的加权（有限）混合。这个滤波器分布的平均值和协方差（假设它们是有限的）可以用bxk表示-1=NXi=1ew（i）k-1bxk-1（i）、（31a）Pk-1=NXi=1ew（i）k-1.主键-1（i）+bxk-1.- bxk-1（i）·T, （31b）其中（A）是AAT的缩写。现在已经观察到了yk，我们想把（26）中的jointposterior传播到时间k，即p（xk）-1，λ1:k-1 | y1:k-1） yk-→ p（xk，λ1:k | y1:k）。（32）这可以通过以下步骤（（1）-（4））实现，如下所述：1）PF预测步骤：从适当选择的方案π（λk |λ（i）1:k生成N个新样本λ（i）kf-1，y1:k）。然后设置λ（i）1:k={λ（i）1:k-1，λ（i）k}，对于i=1，··，N，表示截至时间k的粒子轨迹。关于方案选择，请参见下面的备注4。DRAFT2）KF预测步骤：自之前的p（xk）开始-1 |λ（i）1:k-1，y1:k-1） 由于这一步是高斯的（由（29）给出），并且噪声也是高斯的，由λ（i）k给出（见（25）），所以动态系统现在是线性高斯的，条件是λ（i）1:k。

因此，使用KF，可以解析地获得预测分布，该分布也是高斯分布，并用p（xk |λ（i）1:k，y1:k表示-1） =Nbxk | k-1（i），Pk | k-1（i）, 其中bxk | k-1（i）=Akbxk-1（i）+Bkuwk（i）（33）Pk | k-1（i）=AkPk-1（i）ATk+BkQk（i）BTk。（34）3）KF更新步骤：假设我们现在有了新的观测yk。现在我们可以更新后验分布p（xk |λ（i）1:k，y1:k），由于CLGM，它也是高斯分布，表示为asp（xk |λ（i）1:k，y1:k）=Nbxk（一），Pk（一）. （35）参数bxk（i），Pk（i）可以使用以下步骤从KF递归获得：bxk（i）=bxk | k-1（i）+Kk（i）yk- Ckbxk | k-1（i）- uek（i）（36）Pk（i）=Pk | k-1（i）- Kk（i）CTkPk | k-1（i）、（37）式中kk（i）=Pk | k-1（i）CTkS-1k（i）（38）Sk（i）=CkPk | k-1（i）CTk+Rk（i）。（39）此外，边际似然也以闭合形式获得，它也是高斯的，并且是给定的YP（yk |λ（i）1:k，y1:k-1） =NuLk（i）∑Lk（i）, （40）式中uLk（i）=Ckbxk | k-1（i）+uek（i）（41a）∑Lk（i）=Sk（i）。（41b）4）PF更新步骤：现在给定观测值和粒子{λ（i）1:k}Ni=1，我们需要更新到后验概率（λ1:k | y1:k）=NXi=1ew（i）kδ（λ1:k）- λ（i）1:k）；ew（i）k≥ 0，（42）排水量，PNI=1ew（i）k=1。可以（使用（22））asw（i）k获得相应的重量ew（i）kc∝p（yk |λ（i）1:k，y1:k-1） p（λ（i）k |λ（i）1:k-1，y1:k-1） π（λ（i）k |·）ew（i）k-1，（43a）ew（i）k=w（i）kPNj=1w（j）k.（43b）密度p（yk |λ（i）1:k，y1:k-1） 已在（40）中获得。现在我们假设p（λ（i）k |λ（i）1:k-1，y1:k-1） 是特定的，我们将在IV-B中进一步讨论。这将完成一个传播周期。为了传播到下一个循环，我们首先对粒子λ（i）1:k（以及相关的bxk（一），Pk（一）) 必要时，在再次执行步骤（1）-（4）之前。备注3：只要我们可以直接从WK和ek生成样本，原则上我们可以针对相同的推理任务，例如使用PF。

然而，当状态空间是高维时，众所周知，PF是无效的。另一方面，在我们的RBPF框架中，PF只针对低维状态（辅助混合向量λk）；有条件地使用KF处理剩余的状态向量。因此，RBPF可以很好地适应这里的维度。B.λkAs-PF的动力学规格基本上是为动态演化模型设计的，我们需要指定适当的传递度p（λk |λ1:k）-1，y1:k-1） ，描述路径空间λ1:k的时间演化。然而，这种传递性通常是未知的，我们唯一掌握的信息是λkis根据（已知）平稳密度分布（例如，p*). 我们注意到，之前在对称α稳定噪声的背景下，[46]，[47]用p（λk |λ1:k）处理λkas未知参数-1，y1:k-1） =p（λk）≡ P*（即λkare依赖于k生成）。在PF上下文中，这是完全任意的，从根本上说是很弱的，因为基于路径空间的递归（例如，参见（43a））需要一个转换内核，它将任何新生成的粒子与其祖先血统联系起来。据我们所知，在这种情况下，可证明的态转移密度在文献中没有得到充分的讨论。在我们的方法中，我们通过一个马尔可夫核来指定转移密度，由此核生成的样本构成一个马尔可夫链，其不变分布为边际分布；因为这个不变分布已知，我们用这个不变分布初始化链。在这样做的过程中，我们遵循[48]，作者使用辅助潜在变量构造了AR（1）过程，由此产生的时间序列是平稳的，具有已知的边缘。此外，如原文章所示，时间序列模型遵循简单的自相关函数，该函数呈指数衰减。

这种性质实际上非常有益，因为它可以限制p（xk |λ1:k，y1:k）在粒子路径空间上的充分统计的依赖性。如[48]所示，基于已知的不变边际，当p（λk）属于可完全整除的卷积闭指数族和密度函数族，其中包括（逆）伽马密度、正态密度和β密度时，可以简单地构造这样的转换核。为了保持描述的简单性，而不是深入细节，我们在下面描述一个例子来说明这个想法。关于更多细节，读者可以参考[48]。1）示例：从II-A4中注意，如果X根据多元方差伽马分布，那么边际p（λ）遵循（13）-（14）给出的逆伽马。假设我们对构造p（λk |λk）感兴趣-1). 这可以使用以下关系式[48]λk=ν/2+Ukλk隐式获得-1Vk，（44）英国~ Ga（α，1）和Vk~ Ga（ν/2+α，1），它们相互独立，也相互独立。我们还有（λk |λk）-1) = u(1 - ρ） +ρλk-1，（45）式中ρ=2αν+2α-2是自相关函数，且u=ν-2是边际分布的平均值。注意E（λk |λk-1） 对于ν>2和λkisρ（τ）=ρτ的自相关函数存在。因此，过渡核可以通过一阶AR模型来描述，该模型具有逆伽马边缘。未知参数α（与UK和Vk相关）在这里通过选择的ρ值来确定。注4：π（λk |λ1:k）给出的最优方案[49]-1，y1:k）=π（λk |λk-1，yk）在实践中通常不可用。一个简单的替代方案是用转换核π（λk |λk）实现自举PF-1） 作为提议。C.算法总结算法1总结了针对具有分层高斯噪声的LDS提出的RBPF的自举实现。D