隐马尔科夫模型（HMM）与线性动态系统（LDS）

在我们通常考虑的问题，都是希望假设数据之间是独立的，但有这样一类数据，他们之间有很明显的时序关系，比如天气预报，文字输入等，当把模型建立成独立时，无疑损失了很多有价值的信息。因此我们需要一些模型将时序关系表征在模型里，得到对数据更加准确地描述。HMM和LDS是其中相当经典的模型，在各自领域都有大量应用。
读者可以去我的github阅读一份将HMM用于NER识别的代码，欢迎star，欢迎fork

HMM

马尔可夫模型的缺陷

处理顺序数据的最简单的⽅式是忽略顺序的性质，将观测看做独⽴同分布。然⽽，这种⽅法⽆法利⽤数据中的顺序模式，例如序列中距离较近的观测之间的相关性。例如，假设我们观测⼀个⼆值变量，这个⼆值变量表⽰某⼀天是否下⾬。给定这个变量的⼀系列观测，我们希望预测下⼀天是否会下⾬。如果我们将所有的数据都看成独⽴同分布的，那么我们能够从数据中得到的唯⼀的信息就是⾬天的相对频率。然⽽，在实际⽣活中，我们知道天⽓经常会呈现出持续若⼲天的趋势。因此，观测到今天是否下⾬对于预测明天是否下⾬会有极⼤的帮助。

为了在概率模型中表⽰这种效果，我们需要放松独⽴同分布的假设。完成这件事的⼀种最简单的⽅式是考虑马尔科夫模型（Markov model）。⾸先我们注意到，不失⼀般性，我们可以使⽤概率的乘积规则来表⽰观测序列的联合概率分布，形式为
$$p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_N\right)=p\left({\mathbf{x}}_1\right)\prod _{n=2}^{N}p\left({\mathbf{x}}_n|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}\right)$$

如果我们现在假设右侧的每个条件概率分布只与最近的⼀次观测有关，⽽独⽴于其他所有之前的观测，那么我们就得到了⼀阶马尔科夫链（first-order Markov chain）。这个模型中，N次观测的序列的联合概率分布为
$$p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_N\right)=p\left({\mathbf{x}}_1\right)\sum _{n=2}^{N}p\left({\mathbf{x}}_n|{\mathbf{x}}_{n-1}\right)$$

虽然这⽐独⽴的模型要⼀般⼀些，但是仍然⾮常受限。对于许多顺序的观测来说，我们预计若⼲个连续观测的数据的趋势会为下⼀次预测提供重要的信息。⼀种让更早的观测产⽣影响的⽅法是使⽤⾼阶的马尔科夫链。如果我们允许预测除了与当前观测有关以外，还与当前观测的前⼀次观测有关，那么我们就得到了⼆阶马尔科夫链。然⽽，这种增长的灵活性是有代价的，因为现在模型中参数的数量要多得多。假设观测是具有K个状态的离散变量，那么⼀阶马尔科夫链中的条件概率分布$p\left({\mathbf{x}}_n|{\mathbf{x}}_{n-1}\right)$由K -1个参数指定，每个参数都对应于$x_{n-1}$的K个状态，因此参数的总数为K(K - 1)。现在假设我们将模型推⼴到M阶马尔科夫链，从⽽联合概率分布由条件概率分布$p\left({\mathbf{x}}_n|{\mathbf{x}}_{n-M},\dots ,{\mathbf{x}}_{n-1}\right)$构建。如果变量是离散变量，且条件概率分布使⽤⼀般的条件概率表的形式表⽰，那么这种模型中参数的数量为$K^M(K-1)$。

隐马尔可夫模型

假设我们希望构造任意阶数的不受马尔科夫假设限制的序列模型，同时能够使⽤较少数量的⾃由参数确定。我们可以引⼊额外的潜在变量来使得更丰富的⼀类模型能够从简单的成分中构建，正如我们在第9章讨论混合概率分布和第12章讨论连续潜在变量模型时所做的那样。对于每个观测$x_n$，我们引⼊⼀个对应的潜在变量$z_n$（类型或维度可能与观测变量不同）。我们现在假设潜在变量构成了马尔科夫链，得到的图结构被称为状态空间模型（state space model），它满⾜下⾯的关键的条件独⽴性质，即给定$z_n$的条件下，$z_{n_1}$和$z_{n+1}$是独⽴的，从而
$$z_{n+1}\perp z_{n-1}|z_n$$
于是可以得到观测变量和潜变量的联合分布：
$$p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_N,{\mathbf{z}}_1,\dots ,{\mathbf{z}}_N\right)=p\left({\mathbf{z}}_1\right)\left[\prod _{n=2}^{N}p\left({\mathbf{z}}_n|{\mathbf{z}}_{n-1}\right)\right]\prod _{n=1}^{N}p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)$$
使⽤d-划分准则，我们看到总存在⼀个路径通过潜在变量连接了任意两个观测变量$x_n$和$x_m$，并且这个路径永远不会被阻隔。因此对于观测变量$x_{n+1}$来说，给定所有之前的观测，条件概率分布$p\left({\mathbf{x}}_{n+1}|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)$不会表现出任何的条件独⽴性，因此我们对xn+1的预测依赖于所有之前的观测。
对于顺序数据来说。如果潜在变量是离散的，那么我们得到了隐马尔科夫模型（hidden Markov model）或者HMM（Elliott et al., 1995）。注意，HMM中的观测变量可以是离散的或者是连续的，并且可以使⽤许多不同的条件概率分布进⾏建模。如果潜在变量和观测变量都是⾼斯变量（结点的条件概率分布对于⽗结点的依赖是线性⾼斯的形式），那么我们就得到了线性动态系统（linear dynamical system）。

潜在变量是离散的服从多项分布的变量$z_n$，描述了那个混合分量⽤于⽣成对应的观测$x_n$。⽐较⽅便的做法是使⽤1-of-K表⽰⽅法。 我们现在让$z_n$的概率分布通过条件概率分布$p\left(z_n|z_{n-1}\right)$对前⼀个潜在变量$z_{n_1}$产⽣依赖。由于潜在变量是K维⼆值变量，因此条件概率分布对应于数字组成的表格，记作A，它的元素被称为转移概率（transition probabilities）。元素为$A_{jk}\equiv p\left(z_{nk}=1|z_{n-1,j}=1\right)$。由于它们是概率值，因此满⾜$0\le A_{jk}\le 1$且
${\sum }_k{A}_{jk}=1$，从⽽矩阵A有K(K -1)个独⽴的参数。
这样，我们可以显式地将条件概率分布写成
$$p\left(z_n|z_{n-1},A\right)=\prod _{k=1}^K\prod _{j=1}^K{A}_{jk}^{z_{n-1,j}{z}_{nk}}$$

初始潜在结点$z_1$很特别，因为它没有⽗结点，因此它的边缘概率分布$p(z_1)$由⼀个概率向
量$\pi$表⽰，元素为${\pi }_k\equiv p\left(z_{1k}=1\right)$，即
$$p\left(z_1|\pi \right)=\prod _{k=1}^K{\pi }_k^{z_{1k}}$$

可以通过定义观测变量的条件概率分布p(xn j zn;ϕ)来确定⼀个概率模型，其中ϕ是控制概率分布的参数集合。这些条件概率被称为发射概率（emission probabilities），可以是⾼斯分布（x是连续变量），也可以是条件概率表格（x是离散变量）。由于$x_n$是观测值，因此对于⼀个给定的ϕ值，概率分布$p\left({\mathbf{x}}_n|{\mathbf{z}}_n,\mathbf{\varphi }\right)$由⼀个K维的向量组成，对应于⼆值向量$z_n$的K个可能的状态。我们可以将发射概率表⽰为

$$p\left({\mathbf{x}}_n|{\mathbf{z}}_n,\mathbf{\varphi }\right)=\prod _{k=1}^{K}p{\left({\mathbf{x}}_n|{\mathbf{\varphi }}_k\right)}^{z_{nk}}$$

我们将注意⼒集中在同质的（homogeneous）模型上，其中所有控制潜在变量的条件概率分布都共享相同的参数$A$，类似地所有发射概率分布都共享相同的参数$\varphi$

马尔可夫模型参数的确立

如果我们观测到⼀个数据集 X = { x 1 , … , ∣ x N } \boldsymbol { X } = \left\{ \boldsymbol { x } _ { 1 } , \ldots , | \boldsymbol { x } _ { N } \right\} X={x1​,…,∣xN​}，那么我们可以使⽤最⼤似然法确定HMM的参数。似然函数通过对联合概率分布中的潜在变量进⾏求和的⽅式得到，即
$$p(\mathbf{X}|\mathbf{\theta })=\sum _{Z}p(\mathbf{X},Z|\mathbf{\theta })$$
很自然地，我们希望用EM算法求参数地极大似然。
$$Q\left(\mathbf{\theta },{\mathbf{\theta }}^{旧}\right)=\sum _{Z}p(Z|\mathbf{X},{\mathbf{\theta }}^{旧})\ln p(\mathbf{X},\mathbf{Z}|\mathbf{\theta })$$
为了方便，我们引入两个符号：（注意这两个都是向量）
$$\begin{array}{rl}\gamma \left({\mathbf{z}}_n\right)& =p\left({\mathbf{z}}_n|\mathbf{X},{\mathbf{\theta }}^{旧}\right)\\ \xi \left({\mathbf{z}}_{n-1},{\mathbf{z}}_n\right)& =p\left({\mathbf{z}}_{n-1},{\mathbf{z}}_n|\mathbf{X},{\mathbf{\theta }}^{旧}\right)\end{array}$$
⽤$z_{nk}$来表⽰$z_{nk}=1$的条件概率，类似地使⽤$\xi \left(z_{n-1,j},z_{nk}\right)$来表⽰后⾯介绍的另⼀个概率。
$$\begin{array}{rl}\gamma \left(z_{nk}\right)& =\mathbb{E}\left[z_{nk}\right]=\sum _{z_n}\gamma (z_n)z_{nk}\\ \xi \left(z_{n-1,j},z_{nk}\right)& =\mathbb{E}\left[z_{n-1,j}{z}_{nk}\right]=\sum _{z_{n-1},z_n}\xi \left(z_{n-1},z_n\right)z_{n-1,j}{z}_{nk}\end{array}$$

带入联合分布的表达式：
$$\begin{array}{rl}Q\left(\mathbf{\theta },{\mathbf{\theta }}^{旧}\right)=& \sum _{k=1}^K\gamma \left(z_{1k}\right)\ln {\pi }_k+\sum _{n=1}^N\sum _{j=1}^K\sum _{k=1}^K\xi \left(z_{n-1,j},z_{nk}\right)\ln A_{jk}\\ & +\sum _{n=1}^N\sum _{k=1}^K\gamma \left(z_{nk}\right)\ln p\left({\mathbf{x}}_n|{\varphi }_k\right)\end{array}$$

所以能高效求$\gamma \left({\mathbf{z}}_n\right)$和$\xi \left({\mathbf{z}}_{n-1},{\mathbf{z}}_n\right)$是E步的核心。
下面来详述如何求$\gamma (z_n)$以及$\xi \left({\mathbf{z}}_{n-1},{\mathbf{z}}_n\right)$
首先，我们很容易通过D-划分知道数据的条件独立性满足以下性质：
$$p(\mathbf{X}|{\mathbf{z}}_n)=p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n|{\mathbf{z}}_n\right)p\left({\mathbf{x}}_{n+1},\dots ,{\mathbf{x}}_N|{\mathbf{z}}_n\right)$$
$$p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}|{\mathbf{x}}_n,{\mathbf{z}}_n\right)=p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}|{\mathbf{z}}_n\right)$$
$$p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}|{\mathbf{z}}_{n-1},{\mathbf{z}}_n\right)=p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}|{\mathbf{z}}_{n-1}\right)$$
$$p\left({\mathbf{x}}_{n+1},\dots ,{\mathbf{x}}_N|{\mathbf{z}}_n,{\mathbf{z}}_{n+1}\right)=p\left({\mathbf{x}}_{n+1},\dots ,{\mathbf{x}}_N|{\mathbf{z}}_{n+1}\right)$$
$$p\left({\mathbf{x}}_{n+2},\dots ,{\mathbf{x}}_N|{\mathbf{z}}_{n+1},{\mathbf{x}}_{n+1}\right)=p\left({\mathbf{x}}_{n+2},\dots ,{\mathbf{x}}_N|{\mathbf{z}}_{n+1}\right)$$
$$p(\mathbf{X}|{\mathbf{z}}_{n-1},{\mathbf{z}}_n)=p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}|{\mathbf{z}}_{n-1}\right)p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)p\left({\mathbf{x}}_{n+1},\dots ,{\mathbf{x}}_N|{\mathbf{z}}_n\right)$$
$$\begin{array}{rl}p\left({\mathbf{x}}_{N+1}|\mathbf{X},{\mathbf{z}}_{N+1}\right)& =p\left({\mathbf{x}}_{N+1}|{\mathbf{z}}_{N+1}\right)\\ p\left({\mathbf{z}}_{N+1}|{\mathbf{z}}_N,\mathbf{X}\right)& =p\left({\mathbf{z}}_{N+1},{\mathbf{z}}_N\right)\end{array}$$
首先计算$\gamma (z_{nk})$:
回顾$\gamma (z_n)$的定义, 我们有$\gamma \left({\mathbf{z}}_n\right)=p\left({\mathbf{z}}_n|\mathbf{X}\right)=\frac{p(\mathbf{X}|{\mathbf{z}}_n)p\left({\mathbf{z}}_n\right)}{p(\mathbf{X})}$
我们利用条件独立性容易得：
$$\gamma \left({\mathbf{z}}_n\right)=\frac{p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n,{\mathbf{z}}_n\right)p\left({\mathbf{x}}_{n+1},\dots ,{\mathbf{x}}_N|{\mathbf{z}}_n\right)}{p(\mathbf{X})}=\frac{\alpha \left({\mathbf{z}}_n\right)\beta \left({\mathbf{z}}_n\right)}{p(\mathbf{X})}$$
所以，我们的核心转而放在求$\alpha (z_n)$以及$\beta (z_n)$上面。
我们很容易发现$\alpha (z_n)$有如下性质：
$$\begin{array}{rl}\alpha \left({\mathbf{z}}_n\right)& =p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n,{\mathbf{z}}_n\right)\\ & =p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n|{\mathbf{z}}_n\right)p\left({\mathbf{z}}_n\right)\\ & =p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}|{\mathbf{z}}_n\right)p\left({\mathbf{z}}_n\right)\\ & =p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1},{\mathbf{z}}_n\right)\\ & =p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)\sum _{{\mathbf{z}}_{n-1}}p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1},{\mathbf{z}}_{n-1},{\mathbf{z}}_n\right)\\ & =p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)\sum _{{\mathbf{z}}_{n-1}}p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1},{\mathbf{z}}_n|{\mathbf{z}}_{n-1}\right)p\left({\mathbf{z}}_{n-1}\right)\\ & =p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)\sum _{{\mathbf{z}}_{n-1}}p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}|{\mathbf{z}}_{n-1}\right)p\left({\mathbf{z}}_n|{\mathbf{z}}_{n-1}\right)p\left({\mathbf{z}}_{n-1}\right)\\ & =p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)\sum _{{\mathbf{z}}_{n-1}}p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1},{\mathbf{z}}_{n-1}\right)p\left({\mathbf{z}}_n|{\mathbf{z}}_{n-1}\right)\end{array}$$
所以$\alpha (z_n)$有如下递推关系：
$$\alpha \left({\mathbf{z}}_n\right)=p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)\sum _{{\mathbf{z}}_{n-1}}\alpha \left({\mathbf{z}}_{n-1}\right)p\left({\mathbf{z}}_n|{\mathbf{z}}_{n-1}\right)$$
同理$\beta (z_n)$也有递推关系：
$$\beta \left({\mathbf{z}}_n\right)=\sum _{{\mathbf{z}}_{n+1}}\beta \left({\mathbf{z}}_{n+1}\right)p\left({\mathbf{x}}_{n+1}|{\mathbf{z}}_{n+1}\right)p\left({\mathbf{z}}_{n+1}|{\mathbf{z}}_n\right)$$
如此，我们就能够高效地求解$\alpha (z_n)$与$\beta (z_n)$,进而求出$\gamma (z_n)$

缩放因子

在我们能够在实际应⽤中使⽤前向后向算法之前， 有⼀件事情必须强调，$\alpha (z_n)$很快就会指数地趋近于零。 因此我们使⽤重新缩放的$\alpha (z_n)$和 $\beta (z_n)$来计算，它们的值保持与单位长度在同⼀个量级上。正如我们将看到的那样，当我们在EM算法中使⽤这些缩放的量时，对应的缩放因⼦会消去。
我们考虑：
$$\hat{\alpha }\left({\mathbf{z}}_n\right)=p\left({\mathbf{z}}_n|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)=\frac{\alpha \left({\mathbf{z}}_n\right)}{p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)}$$
令：
$$c_n=p\left({\mathbf{x}}_n|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n-1}\right)$$
根据乘积规则，我们有:
$$p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)=\prod _{m=1}^n{c}_m$$
因此：
$$\alpha \left({\mathbf{z}}_n\right)=p\left({\mathbf{z}}_n|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)=\left(\prod _{m=1}^n{c}_m\right)\hat{\alpha }\left({\mathbf{z}}_n\right)$$
然后我们可以将$\alpha$的递归⽅程（13.36）转化为$\hat{\alpha }$的递归⽅程，形式为
$$c_n\hat{\alpha }\left({\mathbf{z}}_n\right)=p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)\sum _{{\mathbf{z}}_{n-1}}\hat{\alpha }\left({\mathbf{z}}_{n-1}\right)p\left({\mathbf{z}}_n|{\mathbf{z}}_{n-1}\right)$$
同理：
$$c_{n+1}\hat{\beta }\left({\mathbf{z}}_n\right)=\sum _{{\mathbf{z}}_{n+1}}\hat{\beta }\left({\mathbf{z}}_{n+1}\right)p\left({\mathbf{x}}_{n+1}|{\mathbf{z}}_{n+1}\right)p\left({\mathbf{z}}_{n+1}|{\mathbf{z}}_n\right)$$

我们看到所求的边缘概率为
$$\gamma \left({\mathbf{z}}_n\right)=\hat{\alpha }\left({\mathbf{z}}_n\right)\hat{\beta }\left({\mathbf{z}}_n\right)$$
$$\xi \left({\mathbf{z}}_{n-1},{\mathbf{z}}_n\right)=c_n^{-1}\hat{\alpha }\left({\mathbf{z}}_{n-1}\right)p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)p\left({\mathbf{z}}_n|{\mathbf{z}}_{n-1}\right)\hat{\beta }\left({\mathbf{z}}_n\right)$$

利用动态规划高效求解极大似然估计

在隐马尔可夫模型的许多应⽤中，潜在变量有许多有意义的直观意义，因此对于给定的观测序列，我们常常感兴趣的是寻找概率最⾼的隐含状态序列。例如，在语⾳识别中，对于⼀个给定的声⾳观测序列，我们可能希望找到概率最⼤的⾳素序列。在实际应⽤中，我们通常感兴趣的是寻找最可能的状态序列（sequence），这可以使⽤最⼤加和算法⾼效地求出，这个算法在隐马尔科夫模型中被称为维特⽐算法。
事实上，我们很容易注意到，最优路径具有如下特性：如果最优路径在时刻n通过$z_n^{\ast }$。那么从这个点到终点的路径必须是最优的。
我们考虑：
$$\omega \left(z_n\right)=\underset{z_1,\dots ,z_{n-1}}{\max }\ln p\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n,{\mathbf{z}}_1,\dots ,{\mathbf{z}}_n\right)$$
其表达含义为在第n步固定为$z_n$的情况下的前n项对数似然
显然最优解满足$P^{\ast }={\max }_{1\le i\le N}\omega (z_{Ni})$,$z_N^{\ast }=\arg {\max }_i\omega (z_{Ni})$
另外我们需要具体求出 { z 1 , . . . , z N } \{z_1,...,z_N\} {z1​,...,zN​}序列：
事实上我们有：假设定义在时刻n的状态为i的路径中概率最大的n-1个节点为${\varphi }_{n-1}(i)$
$${\varphi }_{n-1}(i)=\arg \underset{j}{\max }[\omega (z_{n-1,j})A_{ji}]$$
于是，我们可以求出每个时刻的潜变量：
$$z_{n-1}^{\ast }={\varphi }_{n-1}(z_n^{\ast })$$

LDS

基本概念

对于顺序数据来说，之前的概率图描述了两个重要的模型。如果潜在变量是离散的，那么我们得到了隐马尔科夫模型（hidden Markov model）或者HMM（Elliott et al., 1995）。注意，HMM中的观测变量可以是离散的或者是连续的，并且可以使⽤许多不同的条件概率分布进⾏建模。如果潜在变量和观测变量都是⾼斯变量（结点的条件概率分布对于⽗结点的依赖是线性⾼斯的形式），那么我们就得到了线性动态系统（linear dynamical system。

为了说明线性动态系统的概念，让我们考虑下⾯这个简单的例⼦，它经常在实际问题中出现。假设我们希望使⽤⼀个有噪声的传感器测量⼀个未知量z的值，传感器返回⼀个观测值x，表⽰z的值加上⼀个零均值的⾼斯噪声。给定⼀个单次的测量，我们关于z的最好的猜测是假设z = x。然⽽，我们可以通过取多次测量然后求平均的⽅法提⾼我们对z的估计效果，因为随机噪声项倾向于彼此抵消。现在，让我们将情况变得更复杂。假设我们希望测量⼀个随着时间变化的量z。我们可以对进⾏常规的测量x，从⽽我们得到了${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N$，我们希望找到对应的$z_1,\dots ,z_N$。如果我们简单地对测量求平均，那么由于随机噪声产⽣的误差会被消去，但是不幸的是我们会仅仅得到⼀个单⼀的平均估计，对z的变化进⾏了平均，从⽽引⼊了⼀种新的误差。

线性动态系统模型

模型假设：
$$\begin{array}{rl}p\left({\mathbf{z}}_n|{\mathbf{z}}_{n-1}\right)& =\mathcal{N}\left({\mathbf{z}}_n|\mathbf{A}{\mathbf{z}}_{n-1},\mathbf{\Gamma }\right)\\ p\left({\mathbf{x}}_n|{\mathbf{z}}_n\right)& =\mathcal{N}\left({\mathbf{x}}_n|{\mathbf{C}\mathbf{z}}_n,\mathbf{\Sigma }\right)\\ p\left(z_1\right)& =\mathcal{N}\left(z_1|{\mu }_0,P_0\right)\end{array}$$
因此，参数为 θ = { A , Γ , C , Σ , μ 0 , P 0 } \boldsymbol { \theta } = \left\{ \boldsymbol { A } , \boldsymbol { \Gamma } , \boldsymbol { C } , \boldsymbol { \Sigma } , \boldsymbol { \mu } _ { 0 } , \boldsymbol { P } _ { 0 } \right\} θ={A,Γ,C,Σ,μ0​,P0​}

线性动态系统的推断：

我们希望考虑$p\left({\mathbf{z}}_n|{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)$
我们不妨将其记作$\hat{\alpha }\left(z_n\right)=\mathcal{N}\left(z_n|{\mathbf{\mu }}_n,{\mathbf{V}}_n\right)$
由于缩放变量$\hat{\alpha }\left(z_n\right)$传播公式完全相同，因此我们有：
$c_n\mathcal{N}\left({\mathbf{z}}_n|{\mathbf{\mu }}_n,{\mathbf{V}}_n\right)=\mathcal{N}\left({\mathbf{x}}_n|\mathbf{C}{\mathbf{z}}_n,\mathbf{\Sigma }\right)\int \mathcal{N}\left({\mathbf{z}}_n|\mathbf{A}{\mathbf{z}}_{n-1},\mathbf{\Gamma }\right)\mathcal{N}\left({\mathbf{z}}_{n-1}|{\mathbf{\mu }}_{n-1},{\mathbf{V}}_{n-1}\right)\mathrm{d}{\mathbf{z}}_{n-1}$
因此，我们有：
$${\mathbf{\mu }}_n=\mathbf{A}{\mathbf{\mu }}_{n-1}+{\mathbf{K}}_n\left({\mathbf{x}}_n-\mathbf{C}\mathbf{A}{\mathbf{\mu }}_{n-1}\right)$$
$${\mathbf{V}}_n=\left(\mathbf{I}-{\mathbf{K}}_n\mathbf{C}\right){\mathbf{P}}_{n-1}$$
后向传播一般考虑$\gamma \left({\mathbf{z}}_n\right)=\hat{\alpha }\left({\mathbf{z}}_n\right)\hat{\beta }\left({\mathbf{z}}_n\right)$

$$\gamma \left({\mathbf{z}}_n\right)=\hat{\alpha }\left({\mathbf{z}}_n\right)\hat{\beta }\left({\mathbf{z}}_n\right)=\mathcal{N}\left({\mathbf{z}}_n|{\hat{\mathbf{\mu }}}_n,{\hat{\mathbf{V}}}_n\right)$$
我们有如下递推关系：
$${\hat{\mathbf{\mu }}}_n={\mathbf{\mu }}_n+{\mathbf{J}}_n\left({\hat{\mathbf{\mu }}}_{n+1}+\mathbf{A}{\mathbf{\mu }}_n\right)$$
$${\hat{\mathbf{V}}}_n={\mathbf{V}}_n+{\mathbf{J}}_n\left({\hat{\mathbf{V}}}_{n+1}-{\mathbf{P}}_n\right){\mathbf{J}}_n^T$$
我们很容易得到$p(Z|X,{\theta }^{旧})$,进而，我们有
$$\mathbb{E}\left[z_n\right]={\hat{\mu }}_n$$
$$\mathbb{E}\left[{\mathbf{z}}_n{\mathbf{z}}_{n-1}^T\right]={\hat{\mathbf{V}}}_n{\mathbf{J}}_{n-1}^T+{\hat{\mathbf{\mu }}}_n{\hat{\mathbf{\mu }}}_{n-1}^T$$
$$\mathbb{E}\left[{\mathbf{z}}_n{\mathbf{z}}_n^T\right]={\hat{\mathbf{V}}}_n+{\hat{\mathbf{\mu }}}_n{\hat{\mathbf{\mu }}}_n^T$$

线性动态模型的参数估计

很自然，我们要用EM算法来求解这个问题：
完全数据的联合似然为：
$$\begin{array}{rl}\ln p(\mathbf{X},Z|\mathbf{\theta })& =\ln p\left({\mathbf{z}}_1|{\mathbf{\mu }}_0,{\mathbf{P}}_0\right)+\sum _{n=2}^N\ln p\left({\mathbf{z}}_n|{\mathbf{z}}_{n-1},\mathbf{A},\mathbf{\Gamma }\right)\\ & +\sum _{n=1}^N\ln p\left({\mathbf{x}}_n|{\mathbf{z}}_n,\mathbf{C},\mathbf{\Sigma }\right)\end{array}$$

$$Q\left(\mathbf{\theta },{\mathbf{\theta }}^{旧}\right)={\mathbb{E}}_{Z|{\mathbf{\theta }}^{旧}}[\ln p(\mathbf{X},Z|\mathbf{\theta })]$$

⾸先考虑参数${\mu }_0$和$P_0$
我们有：
$$Q\left(\mathbf{\theta },{\mathbf{\theta }}^{旧}\right)=-\frac{1}{2}\ln \left|{\mathbf{P}}_0\right|-{\mathbb{E}}_{Z|{\mathbf{\theta }}^{旧}}\left[\frac{1}{2}{\left({\mathbf{z}}_1-{\mathbf{\mu }}_0\right)}^T{\mathbf{P}}_0^{-1}\left({\mathbf{z}}_1-{\mathbf{\mu }}_0\right)\right]+常数$$
$${\mu }_0^{新}=\mathbb{E}\left[z_1\right]$$
$$V_0^{新}=\mathbb{E}\left[{\mathbf{z}}_1{\mathbf{z}}_1^T\right]-\mathbb{E}\left[{\mathbf{z}}_1\right]\mathbb{E}\left[{\mathbf{z}}_1^T\right]$$

类似地可以得到其他参数 θ = { A , Γ , C , Σ } \boldsymbol { \theta } = \left\{ \boldsymbol { A } , \boldsymbol { \Gamma } , \boldsymbol { C } , \boldsymbol { \Sigma }\right\} θ={A,Γ,C,Σ} 的估计

$$\begin{array}{rl}Q\left(\mathbf{\theta },{\mathbf{\theta }}^{旧}\right)& =-\frac{N-1}{2}\ln |\mathbf{\Gamma }|-{\mathbb{E}}_{Z|{\mathbf{\theta }}^{旧}}\left[\frac{1}{2}\sum _{n=2}^N{\left({\mathbf{z}}_n-\mathbf{A}{\mathbf{z}}_{n-1}\right)}^T{\mathbf{\Gamma }}^{-1}\left({\mathbf{z}}_n-\mathbf{A}{\mathbf{z}}_{n-1}\right)\right]\end{array}+常数$$
我们可以求出$A$和$\Gamma$的似然解：
$$A^{新}=\left(\sum _{n=2}^N\mathbb{E}\left[z_n{z}_{n-1}^T\right]\right){\left(\sum _{n=2}^N\mathbb{E}\left[z_{n-1}{z}_{n-1}^T\right]\right)}^{-1}$$
$${\Gamma }^{新}=\frac{1}{N-1}\sum _{n=2}^N\{\mathbb{E}\left[{\mathbf{z}}_n{\mathbf{z}}_n^T\right]-A^{新}\mathbb{E}\left[{\mathbf{z}}_n{\mathbf{z}}_n^T\right]-\mathbb{E}\left[{\mathbf{z}}_n{\mathbf{z}}_n^T\right](A^{新}{)}^T-A^{新}\mathbb{E}\left[{\mathbf{z}}_n{\mathbf{z}}_n^T\right](A^{新}{)}^T$$
以及：
$$\begin{array}{rl}Q\left(\mathbf{\theta },{\mathbf{\theta }}^{\Vert \mathrm{H}}\right)=& -\frac{N}{2}\ln |\mathbf{\Sigma }|-{\mathbb{E}}_{Z|{\mathbf{\theta }}^{旧}}\left[\frac{1}{2}\sum _{n=1}^N{\left({\mathbf{x}}_n-\mathbf{C}{\mathbf{z}}_n\right)}^T{\mathbf{\Sigma }}^{-1}\left({\mathbf{x}}_n-\mathbf{C}{\mathbf{z}}_n\right)\right]\end{array}$$
我们可以求出$C$和$\Sigma$的近似解：
$$C^{新}=\left(\sum _{n=1}^N\mathbb{E}[z_n]\left[z_n^T\right]\right){\left(\sum _{n=1}^N\mathbb{E}\left[z_n{z}_n^T\right]\right)}^{-1}$$
$${\Sigma }^{新}=\frac{1}{N}\sum _{n=1}^N\{\mathbb{E}\left[{\mathbf{x}}_n{\mathbf{x}}_n^T\right]-C^{新}\mathbb{E}\left[{\mathbf{x}}_n{\mathbf{z}}_n^T\right]-\mathbb{E}{\mathbf{x}}_n\left[{\mathbf{z}}_n^T\right](C^{新}{)}^T-C^{新}\mathbb{E}\left[{\mathbf{z}}_n{\mathbf{z}}_n^T\right](C^{新}{)}^T$$