《计算智能》（二）概率

简介

《计算智能》的第二部分——概率，因为各种原因，得分成两部分……

不确定知识与推理

基本概率符号

概率理论中变量被称为随机变量，变量的名字以大写字母开头。每个随机变量有一个定义域，当一个随机变量为布尔变量时，其小写字母表示的是True，例如 $x$ 表示 $X=true$。

当概率 $P$ 表示的是随机变量的分布时，它结果是由一些数组成的向量。

符号 $P$ 也被用于条件分布$P(X|Y)$ 给出每个可能的 $i$、$j$ 组合下的值 $P(X=x_i|Y=y_j)$。

对于连续变量，则 $P$ 为概率密度函数。

除了单个变量的分布外，还需要符号表示多个变量的分布，即为联合概率分布。

使用完全联合分布进行推理

柯尔莫哥洛夫公理：

$$P(a\vee b)=P(a)+P(b)-P(a\wedge b)$$

边缘化规则（求和消元）：

$$P(Y)=\sum _{z\in Z}P(Y,z)$$

条件化规则：

$$P(Y)=\sum _{z}P(Y|z)P(z)$$

通用推理过程：

• 询问变量 $X$，即查询为 $P(X|e)$

• 证据变量 $E$，即观测变量，其观察值为 $e$

• 隐变量，即其余的未观测变量 $Y$

$$P(X|e)=\alpha P(X,e)=\alpha \sum _{y}P(X,e,y)$$

贝叶斯规则：

$$P(Y|X)=\frac{P(X|Y)P(Y)}{P(X)}$$

带证据变量 $e$ 的贝叶斯规则：

$$P(Y|X,e)=\frac{P(X|Y,e)P(Y|e)}{P(X|e)}$$

归一化的贝叶斯规则：

$$P(Y|X)=\alpha P(X|Y)P(Y)$$

给定第三个随机变量 $Z$ 之后，两个随机变量 $X$ 和 $Y$ 的条件独立性的一般定义是：

$$\begin{gathered} P(X,Y|Z)=P(X|Z)P(Y|Z) \\ P(Z|X,Y)=\alpha P(X|Z)P(Y|Z)P(Z) \end{gathered}$$

概率推理

贝叶斯网络

利用全联合概率分布计算不确定知识的局限性：

1. 变量数目增多时，域的规模会增大到不可操作的程度。
2. 为每个原子事件指定概率是很不自然的，而且可能非常困难。

利用贝叶斯网络以一种自然和有效的方式来捕捉不确定知识。该网络是一个用于条件独立断言的简单的图模型，对全联合分布有紧凑的描述形式。

贝叶斯网络，又名：信度网（belief network）、概率网络（probability network）、因果网络（causal network）、知识图（knowledge map）。

定义：贝叶斯网络是一个有向图，其中每个节点都标注了定量概率信息。

• 一个随机变量集组成网络节点。变量可以是离散的或连续的。
• 一个连接点对的有向边或箭头集合。若存在从节点 $X$ 指向节点 $Y$ 的有向边，则称 $X$ 是 $Y$ 的父节点。
• 每个节点 $X_i$ 都有一个条件概率分布：$P(X_i|\text{Parents}(X_i))$
• 一个有向无环图（DAG）

在简单的例子中，条件概率分布被表示为条件概率表（CPT）。CPT中的每一行包含了每个节点值
对于一个条件事件的条件概率。

贝叶斯网络语义

两种方式理解贝叶斯网络语义：

1. 将贝叶斯网络视为对联合分布的表示
2. 将其视为对条件依赖语句集合的编码

联合概率分布通过局部条件概率分布定义：

$$P(X_1,...,X_n)=\prod _{i=1}P(X_i|\text{Parents}(X_i))$$

紧致性与节点排序

贝叶斯网络能够对域进行一种完备而冗余的表示，比全联合概率分布紧凑得多。

• 贝叶斯网络是紧致性的
• 局部结构化（locally structure，sparse）
• 每个组成部分都只与数量有限的其他部分发生直接的相互作用，而不考虑组成部分的总数量。

局部结构的复杂度关于变量 $n$ 是线性的，而不是指数增长的。

• 对布尔变量 $X_i$，有 $k$ 个布尔变量的父节点，CPT中就有 $2^k$个数据。
• 若每个变量的父节点不超过 $k$ 个，则该网络需要的数据是 $O(n·2^k)$
• 即关于 $n$ 是线性的，而全概率分布的数据量是 $O(2^n)$
• 对前面举例的网络,其数据量为 $1+1+4+2+2=10$（vs. $25-1=31$）
• 例如，一个有30个节点的网络，每个节点有5个父节点，那么，该网络需要960个数据，而全联合概率需要10亿个。

构造贝叶斯网络

1. 选择已排序的节点 $X_1,\dots ,X_n$

2. For i=1 to n

   1. 把 $X_i$ 添加到网络

   2. 从 $X_1,\dots ,X_{i-1}$选择父节点

      $$P(X_i|\text{Parents}(X_i))=P(X_i|X_1,...,X_{i-1})$$

      对父节点的选择能保证：

      $$P(X_1,...,X_n)=\prod _{i=1}P(X_i|X_1,...,X_{i-1})=\prod _{i=1}P(X_i|\text{Parents}(X_i))$$

在构造一个局部结构化的域中，

• 不仅要求每个变量只受到少数几个其他变量的直接影响
• 还要求网络的拓扑结构确实反映了合适的父节点集对该变量的那些直接影响

施加直接影响者”将不得不先添加到网络中，因此添加节点的正确次序是首先添加根本原因节点，然后加入受它们直接影响的变量。

贝叶斯网络中的条件独立关系

1. 给定父节点，一个节点与它的非后代节点是条件独立的。Johncall与Burglary和Earthquake是条件独立的。
2. 给定一个节点的父节点、子节点以及子节点的父节点——即，给定它的Markov覆盖，该节点与网络中的所有其它节点都是条件独立的。 Burglary与Johncalls及Marycalls是独立的。

条件分布的有效表达

不确定关系

• CPT随着父节点个数的增加呈指数增长。
• 在父节点或子节点是连续情况下， CPT规模无限增长 。
• 父节点与子节点之间的关系完全任意

解决办法：采用符合某种标准模式的规范分布。 在该情况下，完整的概率分布表能够通过确定所使用的模式，并提供少数几个参数来指定。

确定性节点的提供。一个确定性节点的取值能够由其父节点的取值完全确定，无不确定性。 这种关系可以是一种逻辑关系，也可以是数值的：例如，父节点是几个经销商销售一种特定型号汽车的价格，子节点是一个喜欢还价的人，最后买这种型号汽车所出的价钱，那么子节点就应该是其全部父节点值的最小值。

不确定关系可以用“噪声” 逻辑关系来刻画。噪声或关系（noisy-OR）是逻辑或关系的推广。噪声或模型考虑到每个父节点引起子节点为真的能力的不确定性—父节点与子节点之间的因果关系有可能被抑制。

采用两个假设：

1. 假设所有可能的原因都列出。如果漏掉一些原因，可以增加一个所谓的遗漏节点来涵盖“各种各样的原因”；
2. 它假设每个父节点的抑制独立于其它父节点的抑制。

例如，Fever（发烧）为真，当且仅当Cold（感冒）、Flu（流感）、或者Malaria（疟疾）为真。这样不确定性是指，父结点与子结点之间的因果关系有可能被抑制，因此病人可能得了感冒却没有发烧的症状。给定上述假设，Fever为假当且仅当所有为真的父结点都被抑制，这种情况的概率等于每个父结点的抑制概率的乘积。假设各抑制概率如下：

$$\begin{gathered} q_{\text{cold}}=P(\mathrm{\neg }fever|cold,\mathrm{\neg }flu,\mathrm{\neg }malaria)=0.6 \\ {q}_{\text{flu}}=P(\mathrm{\neg }fever|\mathrm{\neg }cold,flu,\mathrm{\neg }malaria)=0.2 \\ {q}_{\text{malaria}}=P(\mathrm{\neg }fever|\mathrm{\neg }cold,\mathrm{\neg }flu,malaria)=0.1 \end{gathered}$$

那么，根据这个信息以及噪声或的假设，可以建立完整的条件概率表。通用规则是：

$$P(x_i|parents(X_i))=1-\prod _{\{j:X_j=true\}}{q}_j$$

其中的乘积是CPT表中这一行的设置为真的父结点之上的乘积。即

这样，对于其中一个变量依赖于 $k$ 个父节点的噪声逻辑关系，可以用 $O(k)$ 而不是 $O(2^k)$ 个参数来描述其完全条件概率表。

包含连续变量的贝叶斯网络

连续变量具有无限个可能取值 ，连续变量具有无限个可能取值。

解决方法：

1. 离散化随机变量，会导致精度损失，非常巨大的概率表
2. 定义标准的概率密度函数族，通过有限个参数进行指定

混合贝叶斯网络

同时包含离散随机变量和连续随机的网络。

父结点为连续变量

• 对于变量 $Cost$，需要指定 $P(Cost|Harvest,Subsidy)$。离散父节点可以通过直接枚举来处理，即分别指定 $P(Cost|Harvest,subsidy)$ 以及 $P(Cost|Harvest,┐subsidy)$。
• 指定价格 $c$ 依赖于 $Harvest$ 的连续值 $h$ 是如何分布的，也就是将价格分布的参数指定为 $h$ 的一个函数。最常见的选择是线性高斯分布。

当离散变量作为连续变量的父节点加入时，网络就定义一个条件高斯分布：给定全部离散变量的任意赋值，所有连续变量的概率分布是一个多元高斯分布（如上）。

子结点为连续变量

概率单位分布：基本的决策过程具有一个硬阈值，但阈值的精确位置受到随机高斯噪声的影响。

贝叶斯网络中的精确推理

推理任务

计算后验边缘概率 $P(X_i|E=e)$

通过枚举进行推理

不在查询变量和证据变量中的其他变量，即隐变量，须求和。

直接枚举：

$$P(B|j,m)=\alpha \sum _e\sum _{a}P(B)P(e)P(a|B,e)P(j|a)P(m|a)$$

这样对于 $n$ 个布尔变量的网络而言，算法的复杂度为 $O(n·2^n)$。

改进：

$$\begin{array}{rl}P(B|j,m)& =\alpha \sum _e\sum _{a}P(B)P(e)P(a|B,e)P(j|a)P(m|a)\\ & =\alpha P(B)\sum _{e}P(e)\sum _{a}P(a|B,e)P(j|a)P(m|a)\\ & =\alpha [<0.001,0.999>\times 0.002\times <0.95,0.29>\times 0.9\times 0.7]\\ & +<0.001,0.999>\times 0.002\times <0.05,0.71>\times 0.05\times 0.01\\ & +<0.001,0.999>\times 0.998\times <0.94,0.001>\times 0.9\times 0.7\\ & +<0.001,0.999>\times 0.998\times <0.06,0.999>\times 0.05\times 0.01\\ & =\alpha <0.001,0.999>(<0.0012,0.004>+<0.5910,0.0011>)\\ & =\alpha <0.001,0.999><0.5922,0.0015>\\ & =\alpha <0.0005922,0.0014985>\\ & \approx <0.283,0.717>\end{array}$$

枚举法效率低下，存在重复计算。对于一个有 $n$ 个布尔变量的网络，该算法的时间复杂度总是 $O(2^n)$。

变量消元算法

从右到左执行求和，存储中间结果（因子）以避免重新计算。

$$\begin{array}{rl}P(B|j,m)& =\alpha \underset{B}{\underbrace{P(B)}}\sum _e\underset{E}{\underbrace{P(e)}}\sum _a\underset{A}{\underbrace{P(a|B,e)}}\underset{J}{\underbrace{P(j|a)}}\underset{M}{\underbrace{P(m|a)}}\\ & =\alpha P(B)\sum _{e}P(e)\sum _{a}P(a|B,e)P(j|a)f_M(a)\\ & =\alpha P(B)\sum _{e}P(e)\sum _{a}P(a|B,e)f_J(a)f_M(a)\\ & =\alpha P(B)\sum _{e}P(e)\sum _a{f}_A(a,b,e)f_J(a)f_M(a)\\ & =\alpha P(B)\sum _{e}P(e)f_{\bar{A}JM}(b,e)(\text{sum out }A)\\ & =\alpha P(B)f_{\bar{E}\bar{A}JM}(b)(\text{sum out }E)\\ & =\alpha f_B(b)\times f_{\bar{E}\bar{A}JM}(b)\end{array}$$

其中，$f_M(A)=\left[\begin{array}{c}P(m|a)\\ P(m|\mathrm{\neg }a)\end{array}\right]=\left[\begin{array}{c}0.7\\ 0.01\end{array}\right]$，$f_J(A)=\left[\begin{array}{c}P(j|a)\\ P(j|\mathrm{\neg }a)\end{array}\right]=\left[\begin{array}{c}0.9\\ 0.05\end{array}\right]$，$f_A(A,B,E)=[P(a|B,E),P(\mathrm{\neg }a|B,E)]=\left[\begin{array}{cc}0.95& 0.05\\ 0.94& 0.060.29& 0.710.001& 0.999\end{array}\right]$。

$$\begin{array}{rl}{f}_{\bar{A}JM}(B,E)& =\sum _a{f}_A(a,B,E)\times f_J(a)\times f_M(a)\\ & =f_A(a,B,E)\times f_J(a)\times f_M(a)+f_A(\mathrm{\neg }a,B,E)\times f_J(\mathrm{\neg }a)\times f_M(\mathrm{\neg }a)\\ & =\left[\begin{array}{c}0.95\\ 0.94\\ 0.29\\ 0.001\end{array}\right]\times 0.7\times 0.9+\left[\begin{array}{c}0.05\\ 0.06\\ 0.71\\ 0.999\end{array}\right]\times 0.01\times 0.05\\ & =\left[\begin{array}{c}0.5985\\ 0.5922\\ 0.1831\\ 0.0011\end{array}\right]\end{array}$$$$\begin{array}{rl}{f}_{\bar{E}\bar{A}JM}(B)& =f_E(e)\times f_{\bar{A}JM}(B,e)+f_E(\mathrm{\neg }e)\times f_{\mathrm{\neg }AJM}(B,\mathrm{\neg }e)\\ & =0.002\times <0.5985,0.1831>+0.998\times <0.5922,0.0011>\\ & =<0.5922126,0.001464>\end{array}$$$$\begin{array}{rl}P(B|j,m)& =\alpha f_B(B)\times f_{\bar{E}\bar{A}JM}(B)\\ & =\alpha <0.001,0.999><0.5922,0.0015>\\ & =\alpha <0.0005922,0.0014985>\\ & \approx <0.283,0.717>\end{array}$$

从因子乘积中对一个变量进行求和消元是直接计算。注意的技巧是，任何不依赖于将被求和消元的变量都可以移到求和符号的外面。

$$\sum _e{f}_E(e)f_A(A,B,e)f_J(A)f_M(A)=f_J(A)f_M(A)\sum _e{f}_E(e)f_A(A,B,e)=f_J(A)f_M(A)f_{\bar{E},A}(A,B)$$

注意直到我们需要将变量从累积乘积中消去前不会进行矩阵乘法。

考虑另一个查询：$P(JohnCalls|Burglary=true)$。求和式为：

$$P(J|b)=\alpha P(b)\sum _{e}P(e)\sum _{a}P(a|b,e)P(J|a)\sum _{m}P(m|a)$$

其中，$\sum _{m}P(m|a)$等于1，即变量 $M$ 和这个查询无关。可以发现，所有既非查询变量的祖先亦非证据变量的祖先的变量都和查询无关，即可以删除。

精确推理的复杂度

对于单连接网络（或者多树）：如果每个节点的父结点数不超过某个常数，则复杂度与网络节点数呈线性关系。

对于多连接网络：在最坏情况下，具有指数级的时间空间复杂度。

所以可以通过某些算法将多连接网络转化成单连接网络，如：

这类算法称为”聚类算法“，也称为”联合树算法“，也称为”团算法“。

贝叶斯网络中的近似推理

精确推理的局限性：对大规模多连通网络中的精确推理是不可操作的。

近似推理基本思想：

1. 从分布 $S$ 中进行 $N$ 个样本的采样
2. 计算一个逼近的后验概率 $\hat{P}$
3. 证明它收敛到真实概率 $P$

直接采样算法

任何采样算法中最基本的要素是根据已知概率分布生成样本。

直接采样的思想是按照拓扑顺序依次对每个变量进行采样。变量值被采样的概率分布依赖于父结点已得到的赋值。

假设顺序为 $[Cloudy,Sprinkler,Rain,WetGrass]$：

1. 从 $P(Cloudy)=<0.5,0.5>$ 中采样 $Cloudy$，假设返回 $true$。
2. 从 $P(Sprinkler|Cloudy=true)=<0.8,0.2>$ 中采样 $Sprinkler$ ，假设返回 $false$。
3. 从 $P(Rain|Cloudy=true)=<0.1,0.9>$中采样 $Rain$，假设返回 $true$。
4. 从 $P(WetGrass|Sprinkler=false,Rain=true)=<0.9,0.1>$ 中采样 $WetGrass$ ，假设返回 $true$。

这样，PRIOR-SAMPLE（先验采样，即直接采样）返回事件 $[true,false,true,true]$。可以发现，生成的样本服从网络所指定的先验联合概率分布。

令 $S_{PS}(x_1,\dots ,x_n)$为由 Prior-Sample 生成的特定事件的概率。由采样过程得到，

$$S_{PS}(x_1,...,x_n)=\prod _{i=1}^{n}P(x_i|\text{Parents}(X_i))=P(x_1,...,x_n)$$

这个事件的采样概率应该是：

$$S_{PS}(true,false,true,true)=0.5\times 0.9\times 0.8\times 0.9=0.324$$

因此，在 $N$ 的大量样本极限下（ $N$ 趋近于无穷大），我们期望有 $32.4\mathrm{\%}$ 的样本是这个事件。

在任何采样方法中，都是通过对实际生成的样本数来计算答案。假设共有 $N$ 个样本，令 $N(x_1,\dots ,x_n)$ 为特定事件 $x_1,\dots ,x_m$ 发生频率，我们期望该频率在极限情况下收敛到根据采样概率得到的它的期望值：

$$\begin{array}{rl}{\text{lim}}_{N\to \mathrm{\infty }}\hat{P}(x_1,...,x_n)& ={\text{lim}}_{N\to \mathrm{\infty }}{N}_{PS}(x_1,...,x_n)/N\\ & =S_{PS}(x_1,...,x_n)\\ & =P(x_1,...,x_n)\end{array}$$

估计概率在大量样本极限下成为精确值，这样的估计被称为一致估计。

拒绝采样算法

该算法是一种给定一个易于采样的分布，为一个难于采样的分布生成采样样本的通用方法。在其最简单的形式中，它可以用于计算条件概率。

计算步骤：

1. 根据网络指定的先验概率分布生成采样样本
2. 拒绝所有与证据不匹配的样本
3. 在剩余样本中对事件 $X=x$ 的出现频繁程度计数从而得到估计概率 $\hat{P}(X=x|e)$

拒绝采样能够产生真实概率的一致估计。

$$\begin{array}{rl}\hat{P}(X|e)& =\alpha N_{PS}(X,e)\\ & =N_{PS}(X,e)/N_{PS}(e)\\ & \approx P(X,e)/P(e)\\ & =P(X|e)\end{array}$$

假设我们希望采样100个样本来估计 $P(Rain|Sprinkler=true)$。在我们所生成的这100个样本中，假设有73个满足 $Sprinkler=false$，因此被拒绝，同时有27个满足 $Sprinkler=true$；这27个中8个满足 $Rain=true$，19个满足 $Rain=false$。因此，

$$P(Rain|Sprinkler=true)\approx \text{NORMALIZE}(<8,19>)=<0.296,0.704>$$

真实的答案是 $<0.3,0.7>$。随着收集的样本的增多，估计值将收敛到真实值。

局限性：拒绝了太多的样本。随着证据变量个数的增多，与证据 $e$ 相一致的样本在所有样本中所占的比例呈指数级下降，所以对于复杂问题这种方法是完全不可用的。

似然加权算法

只生成与证据 $e$ 一致的事件，避免拒绝采样算法的低效率。

算法描述：该算法固定证据变量 $E$ 的值，只对证据以外的其余变量 $X$ 和 $Y$ 进行采样。 这保证了采样样本与证据一致。在对查询变量的分布进行计数之前，把根据证据得到的事件的似然作为每个事件的权值，该权值通过每个证据变量在给定其父节点取值下的条件概率的乘积进行度量。

求解查询 $P(Rain|Cloudy=true,WetGrass=true)$，假设变量的拓扑顺序是 $Cloudy$，$Sprinkler$、$Rain$、$WetGrass$。过程是这样的：首先，将权值 $w$ 设为1.0。然后生成一个事件：

1. $Cloudy$是一个证据变量，其职位 $true$。因此我们设置$$w\leftarrow w\times P(Cloudy=true)=0.5$$

1. $Sprinkler$ 不是一个证据变量，因此从 $P(Sprinkler|Cloudy=true)=<0.1,0.9>$ 中采样；假设返回 $false$。

2. 类似地，从 $P(Rain|Cloudy=true)=<0.8,0.2>$ 中采样；假设返回 $true$。

3. $WetGrass$ 是一个证据变量，其值为 $true$。因此我们设置

   $$w\leftarrow w\times P(WetGrass=true|Sprinkler=false,Rain=true)=0.45$$

这里似然采样返回权值为0.45的事件 $[true,false,true,true]$，它将被计入 $Rain=true$ 中去。

上述是书上所写的，这里写下自己的理解：似然采样可以看成是，固定住某些变量，即证据变量的直接采样。目标是查询的变量概率，从而对每个事件赋予一定的权值，即概率，说成权值反倒一些不好理解了。按上述所说的，事件 $[t,f,t,t]$ 的权值为0.45，另外还能采样到另外三种事件（总共四种，因为固定了两个证据变量）：事件 $[t,f,f,t]$ 的权值为0、事件 $[t,t,f,t]$ 的权值为 0.45、 事件 $[t,t,t,t]$ 的权值为0.495。查询为，

$$P(Rain|Cloudy=true,WetGrass=true)=\alpha <0.45+0,0.45+0.495>\approx <0.323,0.677>$$

似然加权分析

令 $S_{WS}$ 为加权采样的采样分布，证据变量的固定值为 $e$，证据以外的其他变量为 $Z$。给定父节点后，对 $Z$ 中的每个变量进行采样：

$$S_{WS}(z,e)=\prod _{i=1}^{l}P(z_i|\text{Parents}(Z_i))$$

$\text{Parents}(Z_i)$ 可能包含隐变量和证据变量。$S_{WS}$给予证据关注：每个 $Z_i$ 的采样值会受到它的祖先节点中证据的影响，另外，$S_{WS}$ 对证据的考虑要少于对真实的后验概率 $P(z|e)$ 的考虑。

似然权值 $w$ 补偿了真实采样分布与期望采样分布之间的差距。对已知事件 $z,e$ 的权值

$$w(z,e)=\prod _{i=1}^{m}P(e_i|\text{Parents}(E_i))$$

加权采样概率为：

$$\begin{array}{rl}{S}_{WS}(z,e)w(z,e)& =\prod _{i=1}^{l}P(z_i|\text{Parents}(Z_i))\prod _{i=1}^{m}P(e_i|\text{Parents}(E_i))\\ & =P(z,e)\end{array}$$$$\begin{array}{rl}\hat{P}(x|e)& =\alpha \sum _y{N}_{WS}(x,y,e)w(x,y,e)\\ & \approx {\alpha }^{\prime }\sum _y{S}_{WS}(x,y,e)w(x,y,e)\\ & ={\alpha }^{\prime }\sum _{y}P(x,y,e)={\alpha }^{\prime }P(x,e)\\ & =P(x|e)\end{array}$$

似然加权回到一致估计，使用了所有的样本，效率比拒绝采样高。

局限性：证据变量的个数增加时，它仍然要承受大幅度的性能下降。因为大多数的样本权值都非常低，导致在加权估计中起主导作用的只是那些所占比例很小的、与证据相符合的似然程度不是非常小的样本。

Markov链仿真推理

MCMC：Markov chain Monte Carlo，马尔可夫链蒙特卡洛。

Gibbs采样，一种特殊形式的MCMC算法，特别适合贝叶斯网络。

Gibbs算法描述：

1. 该算法总是通过对前一事件进行随机改变而生成每个事件样本

2. 可以认为网络处于为每个变量指定了值的一个特定的当前状态，而下一个状态则通过非证据变量 $X_i$ 进行采样来产生，其取决于 $X_i$ 的Markov覆盖中的变量的当前值

3. 在状态空间中的随机走动，但保持证据变量的值固定不变。

具体来说，查询 $P(Rain|Sprinkler=T,WetGrass=T)$。证据变量固定为他们的观察值，而隐变量 $Cloudy$ 和 $Rain$ 则随机的初始化：$[T,T,F,T]$。反复执行下面的步骤（对于隐变量和询问变量，即非证据变量，可以以任意顺序）：

1. 对节点 $Cloudy$ 采样，给定它的 $MarkovBlanket$（马尔可夫覆盖）当前值，即根据 $P(Cloudy|Sprinkler=T,Rain=F)$ 采样（父节点、子节点和子节点的父节点）。假设得到采样结果为$Cloudy=false$，即新的状态为 $[F,T,F,T]$。
2. 对节点 $Rain$ 采样，给定它的 $MarkovBlanket$ 当前值，即根据 $P(Rain|Cloudy=F,Sprinkler=T,WetGrass=T)$ 进行采样。假设采样结果为 $Rain=T$。新的当前状态是$[F,T,T,T]$。

这个过程中所访问的每一个状态都是一个样本，能对查询变量 $Rain$ 的估计做贡献。如果该过程访问了20个 $Rain$ 为真的状态和60个 $Rain$ 为假的状态，则所求查询的解为 $NORMALIZE(<20,60>)=<0.25,0.75>$。

MCMC工作机理

基本观点：采样过程最终会进入一种动态平衡。

该特性来自于特定的转移概率，即过程从一种状态转移到另一种状态的概率，通过被采样变量在给定Markov覆盖下的条件概率分布定义。令 $q(x\to x^{\prime })$ 为过程从状态 $x$ 转移到状态 $x^{\prime }$ 的概率，该转移概率定义了状态空间上的Markov链，${\pi }_t(x)$ 为时刻 $t$ 处于状态 $x$ 的概率。${\pi }_{t+1}(x^{\prime })$ 表示在时刻 $t+1$ 处于状态 $x^{\prime }$ 的概率。

给定 ${\pi }_t(x)$，我们对算法可能于时刻 $t$ 到达的所有状态，通过对处于该状态的概率与从该状态转移到状态 $x^{\prime }$ 的概率的乘积求和来计算 ${\pi }_{t+1}(x^{\prime })$：

$${\pi }_{t+1}(x^{\prime })=\sum _x{\pi }_t(x)q(x\to x^{\prime })$$

当 ${\pi }_t(x)={\pi }_{t+1}(x^{\prime })$ 时，Markov链达到了稳态分布，记为 $\pi$，即

$$\pi (x^{\prime })=\sum _x\pi (x)q(x\to x^{\prime })$$

对稳态分布的理解：从每个状态的期望“流出”等于来自于所有状态的期望“流入”。一个明显满足该关系的方式是任何两个状态之间沿两个方向的期望流量相等，称为细致平衡：

$$\pi (x^{\prime })q(x^{\prime }\to x)=\pi (x)q(x\to x^{\prime })$$

细致平衡：

$$\begin{array}{rl}\sum _x\pi (x)q(x\to x^{\prime })& =\sum _x\pi (x^{\prime })q(x^{\prime }\to x)\\ & =\pi (x^{\prime })\sum _{x}q(x^{\prime }\to x)\\ & =\pi (x^{\prime })\end{array}$$

给定Markov覆盖的一个变量独立于其它变量。

$$P(x_i^{\prime }|{\overline{x}}_i,e)=P(x_i^{\prime }|mb(X_i))$$

给定一个变量的概率正比于给定父节点的变量概率与给定各自父节点的每个子节点条件概率的乘积。

$$P(x_i^{\prime }|mb(X_i))=\alpha P(x_i^{\prime }|\text{parents}(X_i))\ast \prod _{Y_j\notin \text{Children}(X_i)}P(y_i|\text{parents}(Y_j))$$

关于时间的概率推理

时间与不确定性

时间片是一个随机变量的集合，其中一部分是可观察的，一部分是不可观察的。

假设所观察到的随机变量属于同一个变量子集：

• $X_t$ 为 $t$ 时刻不可观测的状态变量，状态变量从时刻 $t=0$ 开始；
• $E_t$ 为 $t$ 时刻可观测的证据变量，证据变量从时刻 $t=1$ 开始。

假设时间是离散的，时间间隔依赖于具体问题，

$$X_{a:b}=X_a,X_{a+1},...,X_b$$

转移模型

Markov过程：当前状态只依赖过去有限的已出现状态。最简单的是一阶Markov过程，其中当前状态只依赖于相邻的前一个状态，而与更早的状态无关，即 $P(X_t|X_{0:t-1})=P(X_t|X_{t-1})$，称为转移模型。

观察模型

也称为传感器模型。$P(E_t|X_{0:t},E_{0:t-1})=P(E_t|X_t)$

初始状态模型

$$P(X_0)$$

由上述三个部分能够确定所有变量上完整的联合概率分布：

$$P(X_{0:t},E_{1:t})=P(X_0)\prod _{i=1}^{t}P(X_i|X_{i-1})P(E_i|X_i)$$

若一阶Markov过程不精确，可通过如下方法弥补：

1. 提高Markov过程阶数
2. 扩大状态变量集合。比如，考虑雨季的历史记录，增加变量 $Seaso{n}_t$，或者 $Temperatur{e}_t$，$Humidit{y}_t$，$Pressur{e}_t$。

时序模型中的推理

1. 滤波：$P(X_t|e_{1:t})$

   滤波的任务是计算信念状态，即给定目前为止的所有证据，计算当前状态的后验概率分布。也称为状态估计。滤波是一个理性Agent为掌握当前状态以便进行理性决策所需要采取的行动。例如，给定目前为止对雨伞携带者的过去的所有观察数据，计算今天下雨的概率。

2. 预测：$P(X_{t+k}|e_{1:t})$ 对于 $k>0$

   预测的任务是给定目前为止的所有证据，计算未来状态的后验分布。例如，给定目前为止对雨伞携带者的过去的所有观察数据，计算今天开始三天以后下雨的概率。

3. 平滑：$P(X_k|e_{1:t})$ 对于 $0\le k<t$

   平滑的任务是给定目前为止的所有数据，计算过去某一状态的后验概率。例如，给定目前为止对雨伞携带者的过去的所有观察数据，计算上星期三下雨的概率。平滑为该状态提供了一个比当时能得到的结果更好的估计，因为它结合了更多的证据。

4. 似然解释：${\text{arg max}}_{x_{1:t}}P(x_{1:t}|e_{1:t})$

   给定观察序列，希望找到最可能生成这些观察结果的状态序列。例如，如果前三天每天都出现雨伞，但第四天没出现，那么最可能的解释是前三天下了雨，而第四天没有下。

除了以上推理任务，还有：

• 学习：如果还不知道转移模型和传感器模型，则可以从观察中学习。和静态贝叶斯网络一样，动态贝叶斯网络的学习可以作为推理的一个副产品而完成。推理为哪些转移确实发生和哪些状态会产生传感器读数提供了估计，而且这些估计可以用于对模型进行更新。更新过的模型又提供新的估计，这个过程迭代至收敛。整个算法是期望最大化算法的一个特例。

注意，学习需要的是平滑，而不是滤波，因为平滑提供了对过程状态更好的估计。通过滤波实现的学习可能不会正确地收敛。

滤波

一个有用的滤波算法需要维持一个当前状态估计并进行更新，而不是每次更新时回到整个感知历史。否则，更新代价会随着时间推移越来越大。也就是说，存在某个函数 $f$ 满足：

$$P(X_{t+1}|e_{1:t+1})=f(e_{t+1},P(X_t|e_{1:t}))$$

该过程称为递归估计。可以看成两部分：

1. 将当前的状态分布由时刻 $t$ 向前投影到时刻 $t+1$；
2. 通过新的证据 $e_{t+1}$ 进行更新。

重排公式：

$$\begin{array}{rl}P(X_{t+1}|e_{1:t+1})& =P(X_{t+1}|e_{1:t},e_{t+1})(\mathrm{证}\mathrm{据}\mathrm{分}\mathrm{解})\\ & =\alpha P(e_{t+1}|X_{t+1},e_{1:t})P(X_{t+1}|e_{1:t})(\mathrm{使}\mathrm{用}\mathrm{贝}\mathrm{叶}\mathrm{斯}\mathrm{规}\mathrm{则})\\ & =\alpha P(e_{t+1}|X_{t+1})P(X_{t+1}|e_{1:t})(\mathrm{根}\mathrm{据}\mathrm{传}\mathrm{感}\mathrm{器}\mathrm{马}\mathrm{尔}\mathrm{可}\mathrm{夫}\mathrm{假}\mathrm{设})\end{array}$$

$P(e_{t+1}|X_{t+1})$可从传感器模型，即观察模型得到。而 $P(X_{t+1}|e_{1:t})$ 可通过将当前状态 $X_t$ 条件化，得到下一个状态的单步预测结果：

$$\begin{array}{rl}P(X_{t+1}|e_{1:t+1})& =\alpha P(e_{t+1}|X_{t+1})\sum _{x_t}P(X_{t+1}|x_t,e_{1:t})P(x_t|e_{1:t})\\ & =\alpha P(e_{t+1}|X_{t+1})\sum _{x_t}P(X_{t+1}|x_t)P(x_t|e_{1:t})(\mathrm{马}\mathrm{尔}\mathrm{可}\mathrm{夫}\mathrm{假}\mathrm{设})\end{array}$$

得到递归公式。可以认为滤波估计 $P(X_t|e_{1:t})$ 是向前传播的”消息“ $f_{1:t}$，在每次转移时得到修正，并根据每个新的观察进行更新：

$$f_{1:t+1}=\alpha \text{FORWARD}(f_{1:t},e_{t+1})$$

常数时间，常数空间（与 $t$ 无关）。

求 $P(R_2|u_{1:2})$：

1. 第0天，$P(R_0)=<0.5,0.5>$

2. 第1天，$P(R_1)=\sum _{r_0}P(r_0)P(R_1|r_0)=0.5\times \left(\begin{array}{c}0.70.3\end{array}\right)+0.5\times \left(\begin{array}{c}0.30.7\end{array}\right)=\left(\begin{array}{c}0.50.5\end{array}\right)$

   ​ $P(R_1|u_1)=\alpha P(u_1|R_1)P(R_1)=\alpha \left(\begin{array}{c}0.90.2\end{array}\right)\times \left(\begin{array}{c}0.50.5\end{array}\right)=\alpha \left(\begin{array}{c}0.450.1\end{array}\right)=\left(\begin{array}{c}0.8180.182\end{array}\right)$

3. 第2天，$P(R_2|u_1)=\sum _{r_1}P(R_2|r_1)P(r_1|u_1)=0.818\times \left(\begin{array}{c}0.70.3\end{array}\right)+0.182\times \left(\begin{array}{c}0.30.7\end{array}\right)=\left(\begin{array}{c}0.6270.373\end{array}\right)$

   ​ $P(R_2|u_1,u_2)=\alpha P(u_2|R_2)P(R_2|u_1)=\alpha \left(\begin{array}{c}0.90.2\end{array}\right)\times \left(\begin{array}{c}0.6270.373\end{array}\right)=\alpha \left(\begin{array}{c}0.56430.0746\end{array}\right)=\left(\begin{array}{c}0.8830.117\end{array}\right)$

预测

预测是没有增加新证据的条件下的滤波。

$$\begin{gathered} P(X_{t+k+1}|e_{1:t})=\sum _{x_{t+k}}P(X_{t+k+1}|x_{t+k})P(x_{t+k}|e_{1:t}) \\ P(X_{t+k}|e_{1:t})=\sum _{x_{t+k-1}}P(X_{t+k}|x_{t+k-1})P(x_{t+k-1}|e_{1:t}) \\ ... \end{gathered}$$

类似上述的滤波操作，只是没有新证据能够修正，就只能一直乘转移矩阵……

似然

除了滤波和预测外，还可以利用一种前向递归的方法对证据序列的似然 $P(e_{1:t})$进行计算。这是希望比较可能产生相同证据序列的不同的时序模型，。在这个递归过程中用到一种似然消息 $l_{1:t}(X_t)=P(X_t,e_{1:t})$。不难证明，这个消息的计算与滤波的计算是相同的：

$$l_{1:t+1}=\text{FORWARD}(l_{1:t},e_{t+1})$$

计算出 $l_{1:t}$ 之后，通过求和消元消去 $X_t$ 得到实际似然值：

$$L_{1:t}=P(e_{1:t})=\sum _{x_t}{l}_{1:t}(x_t)$$

不太理解，自己的理解是滤波算的是当前证据变量与当前状态变量的条件概率，而似然算的是当前证据变量与当前状态变量的联合概率，在最后一步求和消元后得到的是，产生这 $t$ 个时刻证据变量的概率，用来比较模型？

平滑

$$\begin{array}{rl}P(X_k|e_{1:t})& =P(X_k|e_{1:k},e_{k+1:t})\\ & =\alpha P(X_k|e_{1:k})P(e_{k+1:t}|X_k,e_{1:t})(\mathrm{使}\mathrm{用}\mathrm{贝}\mathrm{叶}\mathrm{斯}\mathrm{规}\mathrm{则})\\ & =\alpha P(X_k|e_{1:k})P(e_{k+1:t}|X_k)(\mathrm{使}\mathrm{用}\mathrm{条}\mathrm{件}\mathrm{独}\mathrm{立}\mathrm{性})\\ & =\alpha f_{1:k}\times b_{k+1:t}\end{array}$$

其中，$b_{k+1:t}=P(e_{k+1:t}|X_k)$，可以通过后向的递归过程来计算：

$$\begin{array}{rl}P(e_{k+1:t}|X_k)& =\sum _{x_{k+1}}P(e_{k+1:t}|X_k,x_{k+1})P(x_{k+1}|x_k)(\mathrm{将}{X}_{k+1}\mathrm{条}\mathrm{件}\mathrm{化})\\ & =\sum _{x_{k+1}}P(e_{k+1:t}|x_{k+1})P(x_{k+1}|X_k)(\mathrm{根}\mathrm{据}\mathrm{条}\mathrm{件}\mathrm{独}\mathrm{立}\mathrm{性})\\ & =\sum _{x_{k+1}}P(e_{k+1},e_{k+2:t}|x_{k+1})P(x_{k+1}|X_k)\\ & =\sum _{x_{k+1}}P(e_{k+1}|x_{k+1})P(e_{k+2:t}|x_{k+1})P(x_{k+1}|X_k)\end{array}$$

求和式的三个因子中，第一个为观察模型，第三个为转移模型，第二个为递归调用，使用后向的”消息“表示有

$$b_{k+1:t}=\text{BACKWARD}(b_{k+2:t},e_{k+1})$$

和前向递归相同，后向递归中每次更新所需要的时间与空间都是常量，因此与 $t$ 无关。该后向阶段的初始值为 $b_{t+1:t}=P(e_{t+1:t}|X_t)=P(1|X_t)$，其中的 $1$ 表示由1组成的向量。因为 $e_{t+1:t}$ 是一个空序列，观察到它的概率等于1。

给定第1天和第2天都观察到雨伞，要计算 $k=1$ 时下雨概率的平滑估计。由：

$$P(R_1|u_1,u_2)=\alpha P(R_1|u_1)P(u_2|R_1)$$

第一项由前向滤波得到为 $<0.818,0.182>$。而第二项后向递归为：

$$\begin{array}{rl}P(u_2|R_1)& =\sum _{r_2}P(u_2|r_2)P(1|r_2)P(r_2|R_1)\\ & =(0.9\times 1\times <0.7,0.3>)+(0.2\times 1\times <0.3,0.7>)\\ & =<0.69,0.41>\end{array}$$

将其代入公式，发现第1天下雨的平滑估计为：

$$P(R_1|u_1,u_2)=\alpha <0.818,0.182>\times <0.69,0.41>\approx <0.883,0.117>$$

如果我们想平滑整个序列，一个显然的方法就是对每个要平滑的时间步运行一次完整的平滑过程。这导致时间复杂度为 $O(t^2)$。用类似上一节的方法记录已计算的结果，能够降低到 $O(t)$。即记录对整个序列进行前向滤波的结果。然后从时刻 $t$ 到 $1$ 运行后向递归，在每个步骤 $k$，根据已经计算出来的后向消息 $b_{k+1:t}$ 和所存储的前向消息 $f_{1:k}$，计算平滑估计。这个算法称为前向-后向算法。

前向-后向算法局限性：

1. 对于状态空间规模庞大而序列很长的应用，算法的空间复杂度可能会过高。
2. 不适合于联机环境下的计算，在联机环境下算法必须在新证据不断追加到序列末尾的同时，为以前的时间片计算平滑估计。

最大似然解释

给定观察序列，希望找到最可能生成这些观察结果的状态序列。

假设 $[true,true,false,true,true]$ 是警卫观察到的前5天的雨伞序列。解释这个序列的最可能的天气序列是什么？

Viterbi算法：所有时间步上的联合概率。

$$\begin{array}{rl}& {\text{max}}_{x_1,...,x_t}P(x_1,...,x_t,X_{t+1}|e_{1:t+1})\\ & =\alpha P(e_{t+1}|X_{t+1}){\text{max}}_{x_t}(P(X_{t+1}|x_t){\text{max}}_{x_1...x_{t-1}}P(x_1,...,x_{t-1},x_t|e_{1:t}))\end{array}$$

与滤波相比：

1. 前向消息 $f_{1:t}=P(X_t|e_{1:t})$ 被如下消息代替，

   $$m_{1:t}={\text{max}}_{x_1,..,x_{t-1}}P(x_1,...,x_{t-1},X_t|e_{1:t})$$

   即到达每个状态 $x_t$ 的最可能路径的概率；

2. $x_t$ 之上的求和被 $x_t$ 之上的极大值所代替。

上述例子解法：

1. 第1天，$m_1={\text{max}}_{x_0}P(x_0,X_1|e_1)$，因为 $u_1=true$，

   $P(R_1)=\sum _{r_0}P(r_0)P(R_1|r_0)=0.5\times \left(\begin{array}{c}0.70.3\end{array}\right)+0.5\times \left(\begin{array}{c}0.30.7\end{array}\right)=\left(\begin{array}{c}0.50.5\end{array}\right)$

   $P(R_1|u_1)=\alpha P(u_1|R_1)P(R_1)=\alpha \left(\begin{array}{c}0.90.2\end{array}\right)\times \left(\begin{array}{c}0.50.5\end{array}\right)=\alpha \left(\begin{array}{c}0.450.1\end{array}\right)=\left(\begin{array}{c}\mathbf{\text{0.8182}}\mathbf{\text{0.1818}}\end{array}\right)$

2. 第2天，根据上述公式计算，即${\text{max}}_{r_1}(r_1,R_2|u_{1:2})=\alpha P(u_2|R_2){\text{max}}_{r_2}(P(R_2|r_1){\text{max}}_{r_1}P(r_1,r_2|u_{1:2}))$，因为 $u_2=true$，

   $$\begin{array}{rl}& r_1=t,\left(\begin{array}{c}0.9\\ 0.2\end{array}\right)\times \left(\begin{array}{c}0.7\\ 0.3\end{array}\right)\times 0.8182=\left(\begin{array}{c}\mathbf{\text{0.5155}}\\ \mathbf{\text{0.0491}}\end{array}\right)\\ & r_1=f,\left(\begin{array}{c}0.9\\ 0.2\end{array}\right)\times \left(\begin{array}{c}0.3\\ 0.7\end{array}\right)\times 0.1818=\left(\begin{array}{c}0.0327\\ 0.0255\end{array}\right)\end{array}$$

1. 第3天，同上，因为 $u_3=false$，$$\begin{array}{rl}& r_2=t,\left(\begin{array}{c}0.1\\ 0.8\end{array}\right)\times \left(\begin{array}{c}0.7\\ 0.3\end{array}\right)\times 0.5155=\left(\begin{array}{c}\mathbf{\text{0.0361}}\\ \mathbf{\text{0.1237}}\end{array}\right)\\ & r_2=f,\left(\begin{array}{c}0.1\\ 0.8\end{array}\right)\times \left(\begin{array}{c}0.3\\ 0.7\end{array}\right)\times 0.1818=\left(\begin{array}{c}0.0055\\ 0.1018\end{array}\right)\end{array}$$

1. 第4天，同上，因为 $u_4=true$，$$\begin{array}{rl}& r_3=t,\left(\begin{array}{c}0.9\\ 0.2\end{array}\right)\times \left(\begin{array}{c}0.7\\ 0.3\end{array}\right)\times 0.0361=\left(\begin{array}{c}0.0227\\ 0.0022\end{array}\right)\\ & r_3=f,\left(\begin{array}{c}0.9\\ 0.2\end{array}\right)\times \left(\begin{array}{c}0.3\\ 0.7\end{array}\right)\times 0.1237=\left(\begin{array}{c}\mathbf{\text{0.0334}}\\ \mathbf{\text{0.0173}}\end{array}\right)\end{array}$$

1. 第5天，同上，因为 $u_5=true$，$$\begin{array}{rl}& r_4=t,\left(\begin{array}{c}0.9\\ 0.2\end{array}\right)\times \left(\begin{array}{c}0.7\\ 0.3\end{array}\right)\times 0.0334=\left(\begin{array}{c}\mathbf{\text{0.0210}}\\ 0.0020\end{array}\right)\\ & r_4=f,\left(\begin{array}{c}0.9\\ 0.2\end{array}\right)\times \left(\begin{array}{c}0.3\\ 0.7\end{array}\right)\times 0.0173=\left(\begin{array}{c}0.0047\\ \mathbf{\text{0.0024}}\end{array}\right)\end{array}$$

隐马尔可夫模型

简称HMM。

简化的矩阵算法

设状态变量 $X_t$ 的值用整数 $1,\dots ,S$ 表示，其实 $S$ 表示可能状态的数目。转移模型 $P(X_t|X_{t-1})$ 成为一个 $S\times S$ 的矩阵 $T$，其中：

$$T_{ij}=P(X_t=j|X_{t-1}=i)$$

例如，雨伞世界的转移矩阵是：

$$T=P(X_t|X_{t-1})=\left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)$$

传感器模型，即观察模型，也类似。$e_t$ 是已知的，为每个状态指定这个状态使 $e_t$ 出现的概率是多少。将这些值放入一个 $S\times S$ 的矩阵 $O_t$ 中，它的第 $i$ 个对角元素是 $P(e_t|X_t=i)$，其他元素是0。例如，$U_1=true$ 时，$O_1=\left(\begin{array}{ccc}0.9& 00& 0.2\end{array}\right)$；$U_3=false$ 时，$O_3=\left(\begin{array}{cc}0.1& 0\\ 0& 0.8\end{array}\right)$。

前向公式变为：

$$f_{1:t+1}=\alpha O_{t+1}{T}^T{f}_{1:t}$$

后向公式变为：

$$b_{k+1:t}=T{O}_{k+1}{b}_{k+2:t}$$

时间复杂度为 $O(S^2·t)$，空间复杂度为 $O(S·t)$。

前向-后向算法的简单变形

使算法能够在常数空间内执行平滑，而与序列长度无关。

将上述前向公式，向后传递，

$$f_{1:t}={\alpha }^{\prime }(T^T{)}^{-1}{O}_{t+1}^{-1}{f}_{1:t+1}$$

修改后的平滑算法首先执行标准的前向过程以计算 $f_{t:t}$ （抛弃所有中间结果），然后对 $b$ 和 $f$ 同时执行后向过程，用它们来计算每一时间步的平滑估计。

不过这个算法有两个显著的限制：

1. 要求转移矩阵必须是可逆的；

2. 要求传感器模型没有零元素，也就是说，在每个状态下每个观察值都是可能的。

有固定延迟的联机平滑

联机平滑存在一种高效的递归算法——即一种时间复杂度与延迟长度无关的算法。假设延迟为 $d$；要对时间片 $t-d$ 进行平滑，其中 $t$ 表示当前时间。根据平滑公式，需要为时间片 $t-d$ 计算，

$$\alpha f_{1:t-d}\times b_{t-d+1:t}$$

然后，当有了新的观察后，需要为时间片 $t-d+1$ 计算，

$$\alpha f_{1:t-d+1}\times b_{t-d+2:t+1}$$

目标是通过增量方式实现这种计算。首先，由滤波公式，由 $f_{1:t-d}$ 计算 $f_{1:t-d+1}$ 。

因为在旧的后向消息 $b_{t-d+1:t}$ 和新的后向消息 $b_{t-d+2:t+1}$ 之间并不存在简单关系。所以考查旧的后向消息 $b_{t-d+1:t}$ 和序列前端的后向消息 $b_{t+1:t}$ 之间的关系。将矩阵形式的后向公式应用 $d$ 次：

$$b_{t-d+1:t}=(\prod _{i=t-d+1}^{t}T{O}_i)b_{t+1:t}=B_{t-d+1:t}1$$

其中矩阵 $B_{t-d+1:t}$ 为 $T$ 和 $O$ 矩阵序列的乘积。 $B$ 类似“变换算子”，它将后来的后向消息变换成早先的后向消息。类似地，

$$b_{t-d+2:t+1}=(\prod _{i=t-d+2}^{t+1}T{O}_i)b_{t+2:t+1}=B_{t-d+2:t+1}1$$

可以发现，新旧矩阵 $B$ 之间的关系：

$$B_{t-d+2:t+1}=O_{t-d+1}^{-1}{T}^{-1}{B}_{t-d+1:t}T{O}_{t+1}$$

这个公式提供了对 $B$ 矩阵的增量式更新，并进而计算新的后向消息 $b_{t-d+1:t+1}$。

Kalman滤波

隐马尔可夫模型针对的是离散变量；卡尔曼滤波针对的是连续变量。

小鸟的飞行可以用每个时间点的6个连续变量来描述，3个变量用于位置 $(X_t,Y_t,Z_t)$，3个变量用于速度 $({\dot{X}}_t,{\dot{Y}}_t,{\dot{Z}}_t)$。之后使用线性高斯分布来表示转移模型和传感器模型。这意味着下一个状态 $X_{t+1}$ 必须是当前状态 $X_t$ 的线性函数，再加上某个高斯噪声。例如，考虑小鸟的 $X$ 坐标，暂时先忽略其他坐标。令观测之间的间隔为 $\mathrm{\Delta }$，并假设在观测间隔里速度不变；那么位置更新由 $X_{t+\mathrm{\Delta }}=X_t+\dot{X}\mathrm{\Delta }$ 给出。增加高斯噪声（解释风向的变化等）后，得到一个线性高斯转移模型：

$$P(X_{t+\mathrm{\Delta }}=x_{t+\mathrm{\Delta }}|X_t=x_t,{\dot{X}}_t={\dot{x}}_t)=N(x_t+{\dot{x}}_t\mathrm{\Delta },{\sigma }^2)(x_{t+\mathrm{\Delta }})$$

更新高斯分布

滤波：

1. 如果当前分布 $P(X_t|e_{1:t})$ 是高斯分布，并且转移模型 $P(X_{t+1}|x_t)$ 是线性高斯的，那么由

   $$P(X_{t+1}|e_{1:t})={\int }_{X_t}P(X_{t+1}|x_t)P(x_t|e_{1:t})d{x}_t$$

   给出的单步预测分布也是高斯分布。自己理解：从离散扩展到连续，即从求和扩展到积分。

2. 如果预测分布 $P(X_{t+1}|e_{1:t})$ 是高斯分布，传感器模型 $P(e_{t+1}|X_{t+1})$ 是线性高斯的，那么条件化新证据后，更新后的分布

   $$P(X_{t+1}|e_{1:t+1})=\alpha P(e_{t+1}|X_{t+1})P(X_{t+1}|e_{1:t})$$

   也是高斯分布。

因此，卡尔曼滤波的FORWARD算子选取一个高斯前向消息，该消息由均值 ${\mu }_t$ 和协方差矩阵 ${\mathrm{\Sigma }}_t$ 确定；并产生一个新的多元高斯前向消息 $f_{1:t+1}$，该消息由均值 ${\mu }_{t+1}$ 和协方差矩阵 ${\mathrm{\Sigma }}_{t+1}$ 确定。

序列七点为高斯先验概率，

$$f_{1:0}=P(X_0)=N({\mu }_0,{\mathrm{\Sigma }}_0)$$

简单的一维实例

（并不简单。）

假设其先验分布为具有方差 ${\sigma }_0^2$ 的高斯分布：

$$P(x_0)=\alpha e^{-\frac{1}{2}(\frac{(x_0-{\mu }_0{)}^2}{{\sigma }_0^2})}$$

转移模型在当前状态中增加了一个具有常数方差 ${\sigma }_x^2$ 的高斯扰动：

$$P(x_{t+1}|x_t)=\alpha e^{-\frac{1}{2}(\frac{(x_{t+1}-x_t{)}^2}{{\sigma }_x^2})}$$

假设传感器模型具有方差为 ${\sigma }_z^2$ 的高斯噪声：

$$P(z_t|x_t)=\alpha e^{-\frac{1}{2}(\frac{(z_t-x_t{)}^2}{{\sigma }_z^2})}$$

现在，已知先验分布 $P(X_0)$，可以使用上述的滤波第一步计算单步预测分布：

$$\begin{array}{rl}P(x_1)& ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}P(x_1|x_0)P(x_0)d{x}_0\\ & =\alpha {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\frac{1}{2}(\frac{(x_1-x_0{)}^2}{{\sigma }_x^2})}{e}^{-\frac{1}{2}(\frac{(x_0-{\mu }_0{)}^2}{{\sigma }_0^2})}d{x}_0\\ & =\alpha {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\frac{1}{2}(\frac{{\sigma }_0^2(x_1-x_0{)}^2+{\sigma }_x^2(x_0-{\mu }_0{)}^2}{{\sigma }_0^2{\sigma }_x^2})}\end{array}$$

通过配方法将任何二次多项式 $a{x}_0^2+b{x}_0+c$ 改写为平方项 $a(x_0-\frac{-b}{2a}{)}^2$ 与独立于 $x_{)}$ 的余项 $c-\frac{b^2}{4a}$ 之和。余项部分可以从积分中移出，于是得到：

$$P(x_1)=\alpha e^{-\frac{1}{2}(c-\frac{b^2}{4a})}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\frac{1}{2}(a(x_0-\frac{-b}{2a}{)}^2)}d{x}_0$$

现在积分部分就是一个全区间上的高斯分布积分，也就是1.只剩下了二次多项式中的余项。注意到这个预想是关于 $x_1$ 的二次多项式；化简后，

$$P(x_1)=\alpha e^{-\frac{1}{2}(\frac{(x_1-{\mu }_0{)}^2}{{\sigma }_0^2+{\sigma }_x^2})}$$

也就是说，单步预测分布是一个具有相同均值 ${\mu }_0$ 的高斯分布，而其方差等于原来方差 ${\sigma }_0^2$ 与转移方差 ${\sigma }_x^2$ 的和。

之后是第二步，即 $z_1$ 条件化，

$$P(x_1|z_1)=\alpha P(z_1|x_1)P(x_1)=\alpha e^{-\frac{1}{2}(\frac{(z_1-x_1{)}^2}{{\sigma }_z^2})}{e}^{-\frac{1}{2}(\frac{(x_1-{\mu }_0{)}^2}{{\sigma }_0^2+{\sigma }_x^2})}$$

再一次合并指数，并配方，得到，

$$P(x_1|z_1)=\alpha e^{-\frac{1}{2}(\frac{(x_1-\frac{({\sigma }_0^2+{\sigma }_x^2)z_1+{\sigma }_z^2{\mu }_0}{{\sigma }_0^2+{\sigma }_x^2+{\sigma }_z^2}{)}^2}{({\sigma }_0^2+{\sigma }_x^2){\sigma }_z^2/({\sigma }_0^2+{\sigma }_x^2+{\sigma }_z^2)})}$$

由此得到了状态变量的一个新的高斯分布。可以发现，新的均值和标准差可以由原来的均值和标准差按照下面的公式计算得到：

$$\begin{gathered} {\mu }_{t+1}=\frac{({\sigma }_t^2+{\sigma }_x^2)z_1+{\sigma }_z^2{\mu }_t}{{\sigma }_t^2+{\sigma }_x^2+{\sigma }_z^2} \\ {\sigma }_{t+1}^2=\frac{({\sigma }_t^2+{\sigma }_x^2){\sigma }_z^2}{{\sigma }_t^2+{\sigma }_x^2+{\sigma }_z^2} \end{gathered}$$

1. 可以把新均值 ${\mu }_{t+1}$ 解释为新观察 $z_{t+1}$ 和旧均值 ${\mu }_t$ 的一个简单的加权平均。如果观察不可靠，那么 ${\sigma }_z^2$ 很大，则更关注旧均值；如果旧均值不可靠（即 ${\sigma }_t^2$ 很大），或者这个过程高度不可预测（${\sigma }_x^2$ 很大），则更关注观察值。
2. 注意到对方差 ${\sigma }_{t+1}^2$ 的更新是独立于观察的。因此可以在事先计算出方差值的序列
3. 方差值的序列很快收敛到一个固定的值，这个值只与 ${\sigma }_x^2$ 和 ${\sigma }_z^2$ 有关，简化后续计算。

一般情况

(想吐。)

完整的多元高斯分布具有如下形式：

$$N(\mu ,\mathrm{\Sigma })(x)=\alpha e^{-\frac{1}{2}((x-\mu ){\mathrm{\Sigma }}^{-1}(x-\mu ))}$$

首先用卡尔曼滤波定义一般的时序模型。转移模型和传感器模型都允许一个附加高斯噪声的线性变换的线性变换，因此有：

$$\begin{gathered} P(x_{t+1}|x_t)=N(F{x}_t,{\mathrm{\Sigma }}_x)(x_{t+1}) \\ P(z_t|x_t)=N(H{x}_t,{\mathrm{\Sigma }}_z)(z_t) \end{gathered}$$

其中 $F$ 和 ${\mathrm{\Sigma }}_x$是描述线性转移模型和转移噪声协方差的矩阵，而 $H$ 和 ${\mathrm{\Sigma }}_z$ 是传感器模型的相应矩阵。现在的均值与协方差的更新公式为：

$$\begin{gathered} {\mu }_{t+1}=F{\mu }_t+K_{t+1}(z_{t+1}-HF{\mu }_t) \\ {\mathrm{\Sigma }}_{t+1}=(I-K_{t+1}H)(F{\mathrm{\Sigma }}_t{F}^T+{\mathrm{\Sigma }}_x) \end{gathered}$$

其中 $K_{t+1}=(F{\mathrm{\Sigma }}_t{F}^T+{\mathrm{\Sigma }}_x)H^T(H(F{\mathrm{\Sigma }}_t{F}^T+{\mathrm{\Sigma }}_x)H^T+{\mathrm{\Sigma }}_z{)}^{-1}$ 被称为卡尔曼增益矩阵。

扩展的卡尔曼滤波器

卡尔曼滤波器所做的假设是相当强的，线性高斯的转移模型和传感器模型。扩展的卡尔曼滤波器试图克服被建模的系统中的非线性。在一个系统中，如果转移模型不能描述为状态向量的矩阵乘法，这个系统就是非线性的。扩展卡尔曼滤波器的工作机理是将在 $x_t={\mu }_t$ 的区域中的状态 $x_t$ （${\mu }_t$ 是当前状态分布的均值）当作局部线性的，在此基础上对系统进行建模。这对于光滑的、表现良好的系统效果非常好，并且允许追踪者保持和更新一个高斯状态分布（是真实后验概率分布的合理近似）。

对于系统“不平滑”或者“表现不良”的情况，标准的解决办法是采用切换的卡尔曼滤波器。

动态贝叶斯网络

雨伞网络以及卡尔曼滤波器网络都是动态贝叶斯网络的一些例子。

总的来说，动态贝叶斯网络中的每个时间片可以具有任意数量的状态变量 $X_t$ 与证据变量 $E_t$。为了简化，假设变量与有向边从一个时间片到另一个时间片是精确复制的，并且假设动态贝叶斯网络表示的是一个一阶马尔可夫过程，所以每个变量的父结点或者在该变量本身所在那个时间片中，或者在与之相邻的上一个时间片中。

精确推理

每次更新所需要的时间空间复杂度是状态变量个数的指数级，最大因子规模 $O(d^{n+1})$，更新代价 $O(d^{n+2})$。

近似推理

1. 使用样本本身作为当前状态分布的近似表示
2. 将采样集合集中集中于状态空间的高概率区域

粒子滤波

粒子滤波算法的工作机理如下：首先，从先验分布 $P(X_0)$ 中采样得到 $N$ 个初始状体样本构成的总体。然后为每个时间步重复下面的更新循环：

1. 对于每个样本，通过转移模型 $P(X_{t+1}|x_t)$，在给定样本的当前状态值 $x_t$ 条件下，对下一个状态值 $x_{t+1}$ 进行采样，使得该样本前向传播。
2. 对于每个样本，用它赋予新证据的似然值 $P(e_{t+1}|x_{t+1})$ 作为权值。
3. 对总体样本进行重新采样已生成一个新的 $N$ 样本总体。每个新样本是从当前的总体中选取的；某个特定样本被选中的概率与其权值成正比。新样本是没有被赋权值的。

例子：

1. 在时刻 $t$ ，8个样本指示 $Rain$（下面），2个样本指示 $\mathrm{\neg }Rain$（不下雨）。通过转移模型对下一个状态进行采样，每个样本都向前传递。在时刻 $t+1$，6个样本指示 $Rain$，4个样本只是 $\mathrm{\neg }Rain$。
2. 在时刻 $t+1$ 观察到 $\mathrm{\neg }Umbrella$。根据每个样本与这个观察的似然程度给样本赋予权值，图中用圆圈的大小表示。
3. 得到一个10样本的新集合，结果有2个样本指示 $Rain$，8个样本指示$\mathrm{\neg }Rain$。

证明：

假设样本总体开始于在时间 $t$ 的前向消息 $f_{1:t}=P(X_t|e_{1:t})$ 的正确表示。用 $N(x_t|e_{1:t})$ 表示处理完观察 $e_{1:t}$ 之后具有状态 $x_t$ 的样本个数，因此对于足够大的 $N$，有

$$N(x_t|e_{1:t})/N=P(x_t|e_{1:t})$$

现在，在给定时刻 $t$ 的样本的条件下，通过在时刻 $t+1$ 对状态变量进行采样而将每个样本向前传播。从每个 $x_t$ 到达 $x_{t+1}$ 的样本个数等于转移概率乘以 $x_t$ 的总量；因此到达状态 $x_{t+1}$ 的总样本数是

$$N(x_{t+1}|e_{1:t})=\sum _{x_t}P(x_{t+1}|x_t)N(x_t|e_{1:t})$$

现在根据每个样本对于时刻 $t+1$ 的证据的似然为其赋权值。处于状态 $x_{t+1}$ 的样本得到权值 $P(e_{t+1}|x_{t+1})$。因此在观察到证据 $e_{t+1}$ 后处于状态 $x_{t+1}$ 的样本总权值是

$$W(x_{t+1}|e_{1:t+1})=P(e_{t+1}|x_{t+1})N(x_{t+1}|e_{1:t})$$

现在考虑重采样步骤。既然每一个样本都以与其权值成正比的概率被复制，重采样后处于状态 $x_{t+1}$ 的样本数与重新采样前状态 $x_{t+1}$ 中的总权值成正比：

$$\begin{array}{rl}N(x_{t+1}|e_{1:t+1})/N& =\alpha W(x_{t+1}|e_{1:t+1})\\ & =\alpha P(e_{t+1}|x_{t+1})N(x_{t+1}|e_{1:t})\\ & =\alpha P(e_{t+1}|x_{[}t+1])\sum _{x_t}P(x_{t+1}|x_t)N(x_t|e_{1:t})\\ & =\alpha NP(e_{t+1}|x_{t+1})\sum _{x_t}P(x_{t+1}|x_t)P(x_t|e_{1:t})\mathrm{根}\mathrm{据}\mathrm{样}\mathrm{本}\mathrm{数}\mathrm{足}\mathrm{够}\mathrm{大}\mathrm{时}\mathrm{近}\mathrm{似}\mathrm{概}\mathrm{率}\\ & ={\alpha }^{\prime }P(e_{t+1}|x_{t+1})\sum _{x_t}P(x_{t+1}|x_t)P(x_t|e_{1:t})\\ & =P(x_{t+1}|e_{1:t+1})\end{array}$$

因此，经过一轮更新循环后的样本总体正确地表示了时刻 $t+1$ 的前向消息。