马尔可夫链中的期望问题

这个问题是我在做 [ZJOI2013] 抛硬币 - 洛谷这道题的时候了解的一个概念。

在网上也只找到了一篇相关的内容：# 马尔可夫链中的期望问题

故在这里来分享一下其中的期望问题。

马尔可夫链

定义：马尔科夫链为状态空间中经过从一个状态到另一个状态的转换的随机过程。该过程要求具备“无记忆性”，也就是说当前状态的概率仅与上一个状态相关。

数学定义：假定我们的序列状态为 …,Xt−2,Xt−1,Xt,Xt+1,Xt+2,……,Xt−2,Xt−1,Xt,Xt+1,Xt+2,…，那么在 Xt+1Xt+1

时刻的状态的条件概率仅依赖于前一刻的状态 XtXt，即：

P(Xt+1|...,Xt−2,Xt−1,Xt)=P(Xt+1|Xt)P(Xt+1|...,Xt−2,Xt−1,Xt)=P(Xt+1|Xt)

既然某一时刻状态转移的概率只依赖于它的前一个状态 ，那么我们只要能求出系统中任意两个状态之间的转换概率，这个马尔科夫链的模型就定了。

概率转移矩阵

通过马尔科夫链的模型转换，我们可以将事件的状态转换成概率矩阵 （又称状态分布矩阵 ），如下例：

其中线上的数字表示是状态转移到另一状态的概率

额，右下角的应该是D，不想改图了……

那么我们可以构造出如下矩阵：

P=⎛⎝⎜⎜⎜⎜⎜⎜ABCDA0010B0000C00.250.50.25D0.50.500⎞⎠⎟⎟⎟⎟⎟⎟P=(ABCDA0000.5B000.250.5C100.50D000.250)

假定我们初始状态为 AA，于是状态矩阵可以定义为

S=⎛⎝⎜⎜⎜1000⎞⎠⎟⎟⎟S=(1000)

那么经过 kk 次转移后的概率分布是什么？

Sk=PkSSk=PkS

例如 k=4k=4 时

S4=⎛⎝⎜⎜⎜0.06250.1250.250.0625⎞⎠⎟⎟⎟S4=(0.06250.1250.250.0625)

通过 numpy 计算：

s = np.array([[1], [0], [0], [0]])
P = np.array([
      [0, 0, 0, 0.5],
      [0, 0, 0.25, 0.5],
      [1, 0, 0.5, 0],
      [0, 0, 0.25, 0]
])
i, n = 0, 4
while i < n:
    s = P.dot(s)
    i = i + 1
print(s)

转移矩阵的修订

如果我们到达了某一个点就结束了整个过程，那么我们需要对状态转移做一定的更改：

呃，这部分对于本文似乎有一点偏题了……请适当略过

我们新增了一个节点，将所有转移到 BB 的状态转移为最终状态，也就是对于 BB 的状态累计求和。

于是新的矩阵成为：

s = np.array([[1], [0], [0], [0], [0]])
P = np.array([
      [0, 0, 0, 0.5, 0],
      [0, 0, 0.25, 0.5, 0],
      [1, 0, 0.5, 0, 0], 
      [0, 0, 0.25, 0, 0], 
      [0, 1, 0, 0, 1]
])

最终 Sk=PkSSk=PkS 就是积累后的概率矩阵。

状态中的期望

我们关注的是到达状态 BB 所需要的平均转移次数，所以我们把随机变量 XX 看作到达目标状态所需要的步数。

根据变量期望的定义可得：

E(X)=∑k=1∞[k⋅P(X=k)]E(X)=∑k=1∞[k⋅P(X=k)]

那么现在的问题就变为了求 P(X=k)P(X=k)。

根据上文，P(X=k)=[B]Sk=[B]PkSP(X=k)=[B]Sk=[B]PkS，也就是取 SkSk 与 BB 相关的那一项。

或许我们可以求得 [B]Sk[B]Sk 的通项（一个多项式？）从而求得期望。

在顶端链接的文章中提到：

事实上，我们可以计算转移矩阵的 nn 次方，即得到各元素关于 nn 的表达式.

由于我没有学过线性代数，所以我就不在本文中展开。

而是采用第二种方法：期望线性方程组

期望线性方程组

显然至少目前对我而言，求通项的方法不可行，所以我们需要利用另外一种方法。

我们仍然求从 AA 转移到 BB 所需要的平均转移次数。

那么我们不妨设 EiEi 表示从状态 ii 转移到 BB 所需要的期望次数，pjipji 表示从状态 jj 从 ii 转移来的概率。

Ei=1+∑jpji⋅EjEi=1+∑jpji⋅Ej

也就是说 EiEi 等于 ii 可以直接转移到的状态的期望乘上转移的概率的和再加上 11（转移的代价，1步）

以上文中例子为例，可以列出：

EAECED=EC+1=0.25EB+0.25ED+1=0.5EA+0.5EB+1EA=EC+1EC=0.25EB+0.25ED+1ED=0.5EA+0.5EB+1

而特殊的，EB=0EB=0，因为已经到达目标了

而可以解出 EA=4EA=4。

也就是说，从 AA 转移到 BB 的期望步数是 44。

方程矩阵化

如果我们把列出来的方程看作矩阵，那么可以得到：

P = np.array([
    [0-1, 0, 1, 0],
    [0, 0-1, 0, 0],
    [0, 0.25, 0.5-1, 0.25],
    [0.5, 0.5, 0, 0-1]
])

似乎是我原本概率转移矩阵有问题？

也就是原概率矩阵的转置矩阵再在对角线上全部减一。

或者可以表示为

P′=PT−InP′=PT−In

而我们最终求解的方程为：

P′E=⎛⎝⎜⎜⎜⎜−1−1⋮−1⎞⎠⎟⎟⎟⎟P′E=(−1−1⋮−1)

利用死去的高斯消元即可。

链接：[六省联考 2017] 分手是祝愿 - 洛谷

在题解区中全是某种奇怪的 DP 设计。所以这里来讲述一下利用本文中的方法求解的方法。

最开始的模型转化我们就省略了……

我们设 EiEi 表示从剩下 ii 个地方没有按下到全部按下的期望步数。

根据题目，我们有：

Ei=inEi−1+n−inEi+1+1Ei=inEi−1+n−inEi+1+1

并且由于最优限制：

i≤k⟺Ei=ii≤k⟺Ei=i

所以我们可以轻易的利用高斯消元求解……但是很明显，这是不现实的……考虑这个矩阵是一个稀疏矩阵，而且每一个方程是只与临近对角线的 3 个元素相关，所以我们可以轻易的将空间优化到 O(3n)O(3n)。

参考至知乎：# 数值线性代数｜LU分解与稀疏矩阵

于是我们可以旋转一下，只存储非0的元素

接下来我们只需要参考 Guass-Jordan 消元法（先从上到下消下三角，再从下到上消上三角）。那么复杂度成功的也被我们降到了 O(n)O(n)。只是常数有点小大……

作者有话说

其实本文讲述的方法存在很大的局限性，尤其的求解的过程，很容易就变成了 O(n3)O(n3)。

不过似乎建模的过程还是蛮简单的？

在求期望的问题中，这不失为一个思考的方向。有些时候，可以通过题目相关的特殊性质进行优化，使得复杂度降低。

可是毒瘤出题人们……还是看每道题的正解吧，或许终将有一天，这个方法可以被一种通法优化到 O(n2)O(n2) 呢？

算法学习笔记(23): 马尔可夫链中的期望问题 - jeefy

马尔可夫链中的期望问题