---

前言

第一節課的后半段其實就是一個馬爾可夫的實際案例教學，我這里在網上找到一個合適案例進行學習，cs234 的課程感覺有點空。

7 markov 實踐

7.1 markov 過程再敘

markov 過程就是一個狀態轉移過程，且該當前狀態只和上一個狀態有關，和歷史無關。
即：$P(s_t|s_{t?1})=P(s_t|s_{t?1},s_{t?2},s_{t?3},...,s_n)$

markov 狀態轉移矩陣：
$$P=\begin{array}{c}\left[\begin{array}{cccc}{p}_{s_1|s_1}& p_{s_2|s_1}& ...& p_{s_n|s_1}\\ p_{s_2|s_1}& p_{s_2|s_2}& ...& p_{s_n|s_2}\\ ...& & & \\ p_{s_1|s_n}& p_{s_2|s_n}& ...& p_{s_n|s_n}\end{array}\right]\end{array}\text{\hspace{1em}?}$$
矩陣中的 第i行第j列表示狀態$s_i$ 到 $s_j$ 的概率$p(s_i|s_j)=p(s_{t+1}=s_j|s_t=s_i)$。稱
$P(s^{\prime }|s)$ 為轉移函數。這里要求 從某個狀態到其他所有狀態的概率和必須為1， 即P矩陣每行概率和為1

如果我們按照滿足markov 性，根據狀態轉移矩陣，得到一個狀態轉移序列
s1->s1->s2->s3->s4->s5->s6 那么就得到了一個markov chain 即馬爾可夫鏈。

7.2 markov 獎勵過程 MRP（markov reward process）

馬爾可夫獎勵過程是 [S, P r, $\gamma$]
S：狀態集合
P：狀態轉移矩陣
r ： reward
$\gamma$：discount factor

為什么要獎勵（reward）？
（1）一個穩定的世界需要反饋，合理的反饋可以讓我們趨于一個穩定。
因此引入獎勵機制。我們將獎勵機制和markov 過程結合，那么有
（2）我們針對不同場景，有不同的回報，因此獎勵機制可以調整我們如何適應變化的環境。

為什么要折扣（discount factor）
（1）一個馬爾可夫過程是有可能出現閉環，如果無限循環下去，那么獎勵就有可能無限累加，要避免這種獎勵因子不斷累加，那么就需要折扣。在想想這句古話：一股做氣，再而衰，三而竭。這不就是折扣因子么。
（2）有時候我們需要近期的效果，那么我們會將長遠利益打一些折扣。相反，我們關注長遠利益時，需要近期利益打折扣

將reward 和 discount factor 結合得到回報（Return）
$G_t=R_t+\gamma ?R_{t+1}+{\gamma }^2?R_{t+2}+...={\sum }_{k=0}^{\infty }{\gamma }^k{R}_{t+k}$

7.3 markov 價值函數與貝爾曼方程

價值（value）：一個狀態的期望回報，即從這個狀態出發的未來累積獎勵的期望值，被稱為這個狀態的價值:

$$\begin{gathered} V(s) \\ =E[R_t+\gamma ?R_{t+1}+{\gamma }^2?R_{t+2}+...|s=s_t] \\ =E[R_t+\gamma (R_{t+1}+\gamma ?R_{t+2}+...|s=s_t)] \\ =E[R_t+\gamma G(s=s_{t+1}|s=s_t)] \\ =E[R_t+\gamma G_{t+1}|s=s_t] \\ =E[R_t|s_t]+E\gamma G_{t+1}|s_t \\ =r(s)+\gamma V_{t+1}|s=s_t \\ =r(s)+\gamma \sum p(s_{t+1}|s_t)V_{t+1}，注：我第一遍推成r(s)+\gamma \sum p(s_{t+1}|s_t)V_t，導致后面直接推不下去了 \\ =r(s)+\gamma \sum p(s^{\prime }|s)V(s^{\prime }),s^{\prime }\in S \\ =貝爾曼方程（BellmanEquation） \end{gathered}$$

于是我們不難得到：
當s’= s1的時候：
$$\begin{array}{c}\left[\begin{array}{c}V(s_1)\end{array}\right]=\left[\begin{array}{c}r(s_1)\end{array}\right]+\gamma \left[\begin{array}{cccc}{p}_{s_1|s_1}& p_{s_2|s_1}& ...& p_{s_n|s_1}\end{array}\right]\left[\begin{array}{c}V(s1)\\ V(s2)\\ ...\\ V(sn)\end{array}\right]\end{array}\text{\hspace{1em}?}$$
當s’=s2的時候：
$$\begin{array}{c}\left[\begin{array}{c}V(s_2)\end{array}\right]=\left[\begin{array}{c}r(s_2)\end{array}\right]+\gamma \left[\begin{array}{cccc}{p}_{s_1|s_2}& p_{s_2|s_2}& ...& p_{s_n|s_2}\end{array}\right]\left[\begin{array}{c}V(s1)\\ V(s2)\\ ...\\ V(sn)\end{array}\right]\end{array}\text{\hspace{1em}?}$$

接下來，我們寫成矩陣形式：

$$\begin{array}{c}\left[\begin{array}{c}V(s1)\\ V(s2)\\ ...\\ V(sn)\end{array}\right]=\left[\begin{array}{c}r(s1)\\ r(s2)\\ ...\\ r(sn)\end{array}\right]+\gamma \left[\begin{array}{cccc}{p}_{s_1|s_1}& p_{s_2|s_1}& ...& p_{s_n|s_1}\\ p_{s_2|s_1}& p_{s_2|s_2}& ...& p_{s_n|s_2}\\ ...& & & \\ p_{s_1|s_n}& p_{s_2|s_n}& ...& p_{s_n|s_n}\end{array}\right]\left[\begin{array}{c}V(s1)\\ V(s2)\\ ...\\ V(sn)\end{array}\right]\end{array}\text{\hspace{1em}?}$$

于是就得到：
$V=R+\gamma PV$
$V?\gamma PV=R$
$(I?\gamma P)V=R$
$V=(I?\gamma P{)}^{?1}R$

按照以往計算經驗，這個矩陣解起來巨麻煩，所以會用 動態規劃(dynamic programming)、 **蒙特卡羅(模擬特-Carlo method)**方法 或 時序差分（temporal difference）

7.4 markov 決策過程MDP（markov decision process）的 狀態價值函數

7.4.1 狀態價值函數

智能體（agent）的策略（Policy）通常用$\pi$表示。策略$\pi (a|s)=P(A_t=a|S_t=s)$是一個函數，表示在s狀態下采取a動作的概率。當一個策略是確定性策略(deterministic policy)的時候，那么智能體在每個狀態只輸出一個確定動作。
當智能體的策略是隨機測策略（stochastic policy）時，那么這個函數輸出的是關于動作的概率分布。

狀態價值函數:
我們用$V^{\pi }(s)$表示在MDP基于測率$\pi$策略得到的價值函數期望：
$V^{\pi }(s)=E_{\pi }[G_t|S_t=s]$，
我這里專門推敲了下：$V^{\pi }$和 V 是一回事，只是為了講名是什么策略，因此加了$\pi$，即乘以一個概率。

7.4.2 狀態價值函數的 貝爾曼期望方程（Bellman Expectation Equation）

根據上面貝爾曼方程算$V(s)$：
$V(s)=r(s)+\gamma \sum p(s^{\prime }|s)V(s^{\prime }),s^{\prime }\in S$
當我們要將在那個策略下時，不難得到:
$$\begin{gathered} V^{\pi }(s)=E_{\pi }[G_t|S_t=s] \\ =\sum \pi (a|s)[r(s)+\gamma \sum p(s^{\prime }|s)V(s^{\prime })],s^{\prime }\in S \end{gathered}$$——因為需要策略$\pi$得概率，因此需要乘以$\pi$

$V^{\pi }(s)=r(s,a)+\gamma {\sum }_{s^{\prime }\in S}p(s^{\prime }|s,a){\sum }_{a^{\prime }\in A}\pi (a|s^{\prime })Q^{\pi }(s^{\prime }，a^{\prime })$

7.5 markov 決策過程MDP（markov decision process）的 動作價值函數

markov 決策過程MDP（markov decision process）—— 動作價值函數

7.5.1 動作價值函數

不同于MRP，在MDP過程中，由于動作的存在，額外定義一個動作價值函數（action-value function）。用$Q^{\pi }$ 表示，在s 狀態下，執行動作a的得到的期望：
$Q^{\pi }(s,a)=E_{\pi }[G_t|S_t=s,A_t=a]$ 。說實話我這里被定義給搞暈了，因此我理解這里就是不需要乘以$\pi (a|s)$

7.5.2 動作價值函數和狀態價值函數

所以得到$V^{\pi }$和$Q^{\pi }$的關系：

（1）$V^{\pi }(s)={\sum }_{a\in A}\pi (a|s)Q^{\pi }(s,a)$
這個式子描述的是：使用策略$\pi$, 狀態 s的價值期望，等于動作價值函數乘以發生動作概率的乘積的總和。這里是動作未發生，需要乘上動作的概率和動作的價值

（2）$Q^{\pi }(s,a)=r(s,a)+\gamma \sum P(s^{\prime }|s,a)V^{\pi }(s^{\prime })$
使用策略$\pi$時，狀態s下采取a動作后的價值期望 等于 當下的獎勵加上 經過$\gamma$衰減之后的所有狀態狀態轉移概率與相應價值的乘積。
這里是動作已經確定，但是狀態不確定，因此乘的是狀態轉移矩陣和狀態

也就是說狀態與狀態之間不再是單純的轉移，還有動作的這個價值反饋加進去。

7.5.3 動作價值函數的貝爾曼期望等式

根據定義：
$Q^{\pi }(s,a)=r(s,a)+\gamma {\sum }_{s^{\prime }\in S}P(s^{\prime }|s,a)V^{\pi }(s^{\prime })$

又因為：
$V^{\pi }(s)={\sum }_{a\in A}\pi (a|s)Q^{\pi }(s,a)$
那么：
$V^{\pi }(s^{\prime })={\sum }_{a\in A}\pi (a|s^{\prime })Q^{\pi }(s^{\prime },a)$
上面的式子可以再變個型，帶入后：
$$\begin{gathered} Q^{\pi }(s,a)=r(s,a)+\gamma {\sum }_{s^{\prime }\in S}P(s^{\prime }|s,a)V^{\pi }(s^{\prime }) \\ =r(s,a)+\gamma {\sum }_{s^{\prime }\in S}P(s^{\prime }|s,a){\sum }_{a\in A}\pi (a|s^{\prime })Q^{\pi }(s^{\prime },a) \end{gathered}$$
——————————————————————
動作價值函數和狀態價值函數的貝爾曼方程很常見，所以我這里推敲了下。
——————————————————————

7.6 動作價值和狀態價值例子

圖中的
（1）虛線表示動作到狀態
（2）圖中的實現表示從當前狀態開始當前動作
（3）紅色的數字是標記狀態獎勵
（4）沒有標記數字的線表示概率為1,如果標記了表示對應概率。

7.6.1 邊緣化（marginalization)，就

這里有個一般的計算MDP的方法，就是將測率的動作邊緣化:
(1) 得到一個沒有動作的mrp。即對于某個狀態，我跟將根據動作策略進行加權，就得到r’(s)是該狀態下的獎勵：
$s^{\prime }={\sum }_{a\in A}\pi (a|s)r(s,a)$

(2)同理，將計算采取動作的概率$\pi$與 將s轉移到s’轉移矩陣進行相乘再累加，就得到一個MRP的從s轉移到s‘的轉移概率。

這樣的做法有如下好處：
簡化問題結構
MDP 涉及狀態和動作兩層結構，分析和求解復雜。而給定策略后，動作選擇就變成確定的概率分布，這時只剩下狀態和狀態之間的轉移（+獎勵）——正好就是 MRP 的結構。
MRP 更容易求解
MRP 沒有動作維度，可以直接用線性代數（比如矩陣解法或貝爾曼方程迭代）來求解狀態值函數 𝑉(𝑠)非常方便。

這個方法是*給定策略下的通用方法。 如果你要做最優策略求解（如值迭代、Q-learning），那就不能只轉成 MRP，因為你需要在每一步決策中“尋找最優動作”，那就是另一套框架了（比如 Bellman Optimality Equation）。

實踐代碼如下：

import numpy as npdef join(str1, str2):return str1 + '-' + str2# markov_chain is a list that save the index of state 
# for example:s1->s1->s2->s3->s4->s5->s6 
# start_index: the start index for markov_chain, not the state index.
def get_return(start_index, markov_chain, rewards, gamma):G_t = 0# >>>>>>> the code is very tricky but effective! <<<<<<<for idx in reversed(range(start_index, len(markov_chain))):# the state index is start at 1 end in 6, so when we use the idx we need to minus 1G_t = gamma * G_t + rewards[markov_chain[idx] - 1]return G_tdef get_value(p_matrix, rewards, gamma):states_amount = len(rewards)rewards = np.array(rewards)rewards = rewards.reshape((-1, 1))value = np.dot(np.linalg.inv(np.eye(states_amount, states_amount) - gamma * p_matrix), rewards)return valuedef mrp_test():np.random.seed(0)p = [[0.8, 0.1, 0.1, 0.0, 0.0, 0.0],[0.0, 0.2, 0.5, 0.3, 0.0, 0.0],[0.0, 0.0, 0.5, 0.5, 0.0, 0.0],[0.1, 0.1, 0.1, 0.1, 0.3, 0.3],[0.2, 0.5, 0.1, 0.0, 0.1, 0.1],[0.0, 0.0, 0.2, 0.3, 0.4, 0.1],]p = np.array(p)# fist to checkt if the sum of each row is 1.0H, W = p.shapefor h in range(H):sum_h = np.sum(p[h,:])# for float compoare we can not suppose it will be exactly 1.0 actually is 0.999999999....if sum_h < 0.9999999 or sum_h > 1.0: print("error in line:" + str(h) + " sum:" + str(sum_h))exit()#          s1, s2, s3, s4, s5, s6rewards = [-1, -2,  0,  1,  2,  4] gamma = 0.7markov_chain = [1, 1, 2, 3, 4, 5, 6, 1]# get the return for agentG_t = get_return(0, markov_chain, rewards, gamma)# get the value for agentV_t = get_value(p, rewards, gamma)print(">>> markov finish!")def mdp_test():states = ["s1", "s2", "s3", "s4", "s5", "s6"]actions = ["hold_s1", "arrival_s2","arrival_s3","arrival_s4", "arrival_s5", "arrival_s6", "stochastic_arrival"]p = {"s1-hold_s1-s1": 1.0,"s1-arrival_s2-s2": 1.0,"s2-arrival_s1-s1": 1.0,"s2-arrival_s3-s3": 1.0,"s3-arrival_s4-s4": 1.0,"s3-arrival_s6-s6": 1.0,"s4-arrival_s5-s5": 1.0,"s5-arrival_s6-s6": 1.0,"s4-stochastic_arrival_s2": 0.2,"s4-stochastic_arrival_s3": 0.3,"s4-stochastic_arrival_s4": 0.5,}rewards = {"s1-hold_s1":   -1,"s1-arrival_s2": 0,"s2-arrival_s1":-1,"s2-arrival_s3":-2,"s3-arrival_s4":-2,"s3-arrival_s6": 7,"s4-arrival_s5": 2,"s5-arrival_s6": 8,"s4-stochastic_arrival_s2": 1,"s4-stochastic_arrival_s3": 1,"s4-stochastic_arrival_s4": 1 }gamma = 0.5mdp = (states, actions, p, gamma)pi_1 = {"s1-hold_s1": 0.5,"s1-arrival_s2": 0.5,"s2-arrival_s1": 0.5,"s2-arrival_s3": 0.5,"s3-arrival_s4": 0.5,"s3-arrival_s5": 0.5,"s4-arrival_s5": 0.5,"s4-stochastic_arrival": 0.5,} pi_2 = {"s1-hold_s1": 0.6,"s1-arrival_s2": 0.4,"s2-arrival_s1": 0.3,"s2-arrival_s3": 0.7,"s3-arrival_s4": 0.5,"s3-arrival_s5": 0.5,"s4-arrival_s5": 0.1,"s4-stochastic_arrival": 0.9,} # 轉化后的MRP的狀態轉移矩陣p_mdp2mrp_pi_1 = [# s1, s2,  s3,  s4,  s5,  s6[1.0 * 0.5, 1.0 * 0.5, 0.0,       0.0,       0.0,       0.0      ],[1.0 * 0.5, 0.0,       1.0 * 0.5, 0.0,       0.0,       0.0      ],[0.0,       0.0,       0.0,       1.0 * 0.5, 0.0,       1.0 * 0.5],[0.0,       0.2 * 0.5, 0.3 * 0.5, 0.5,       1.0 * 0.5, 0.0      ],[0.0,       0.0,       0.0,       0.0,       0.0,       1.0 * 0.5],[0.0,       0.0,       0.0,       0.0,       0.0,       1.0      ],]p_mdp2mrp_pi_1 = np.array(p_mdp2mrp_pi_1)R_mdp2mrp_pi_1 = [-1 * 0.5, -1 * 0.5 + -2 * 0.5,7 * 0.5 + -2 * 0.5,2 * 0.5 +  1 * 0.5, 8 * 0.5,0]# get the  mrp base on mdpv = get_value(p_mdp2mrp_pi_1, R_mdp2mrp_pi_1, gamma)print("mdp v=" + str(v))print("mdp finish!")if __name__ == "__main__":#mrp_test()mdp_test()

運行結果：