mirror of
https://github.com/microsoft/ai-edu.git
synced 2026-08-28 17:41:25 +08:00
Xiaowu/20020407 (#722)
* more * add new * upd * jjk * add drive * Update 安全驾驶问题.md * upadte * fix * add 42 * update * update * update * update * update * rename * asd * Update 贝尔曼方程.md * upadte * update * Update 贝尔曼方程.md * update
This commit is contained in:
xiaowuhu
@@ -85,40 +85,36 @@ RMSE(a,y)=\sqrt{[(0.1-0)^2 + (1-1)^2 + (1.8-2)^2)]/3}=0.129
|
||||
RMSE(b,y)=\sqrt{[(0.1-0)^2 + (1.1-1)^2 + (1.8-2)^2)]/3}=0.141
|
||||
$$
|
||||
|
||||
MC-1
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
G_1 &= R_1
|
||||
\\
|
||||
G_2 &= G_1 + \gamma R_2=R_1+\gamma R_2
|
||||
\\
|
||||
G_3 &= G_2 + \gamma^2 R_3 = R_1 + \gamma R_2 + \gamma^2 R_3
|
||||
\\
|
||||
G_4 &= G_3 + \gamma^3 R_4 = R_1 + \gamma R_2 + \gamma^2 R_3 + \gamma^3 R_4
|
||||
G_{[1]} &= R_1
|
||||
\\\\
|
||||
G_{[2]} &= G_{[1]} + \gamma R_2=R_1+\gamma R_2
|
||||
\\\\
|
||||
G_{[3]} &= G_{[2]} + \gamma^2 R_3 = R_1 + \gamma R_2 + \gamma^2 R_3
|
||||
\\\\
|
||||
G_{[4]} &= G_{[3]} + \gamma^3 R_4 = R_1 + \gamma R_2 + \gamma^2 R_3 + \gamma^3 R_4
|
||||
\\\\
|
||||
G_{[5]} &= G_{[4]} + \gamma^4 R_T = R_1 + \gamma R_2 + \gamma^2 R_3 + \gamma^3 R_4+ \gamma^4 R_T, & V[S_s] += G_{[5]}
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
G[4] &= R_4
|
||||
\\
|
||||
G[3] &= \gamma G[4]+R_3 = R_3 + \gamma R_4
|
||||
\\
|
||||
G[2] &= \gamma G[3]+R_2 = R_2 + \gamma R_3 + \gamma^2 R_4
|
||||
\\
|
||||
G[1] &= \gamma G[2]+R_1 = R_1 + \gamma R_2 + \gamma^2 R_3 + \gamma^3 R_4
|
||||
\end{aligned}
|
||||
$$
|
||||
MC-2
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
G_1 &= R_4, & V[S_{9}]+=G_1
|
||||
\\
|
||||
G_2 &= \gamma G_1 + R_3 = R_3 + \gamma R_4, & V[S_5]+=G_2
|
||||
\\
|
||||
G_3 &= \gamma G_2 + R_2 = R_2 + \gamma R_3 + \gamma^2 R_4, & V[S_4]+=G_3
|
||||
\\
|
||||
G_4 &= \gamma G_3 + R_1 = R_1 + \gamma R_2 + \gamma^2 R_3 + \gamma^3 R_4 , & V[S_0] += G_4
|
||||
G_{[5]} &= \gamma G_{[4]} + R_1 = R_1 + \gamma R_2 + \gamma^2 R_3 + \gamma^3 R_4+\gamma^4 R_T, & V[S_S] += G_{[5]}
|
||||
\\\\
|
||||
G_{[4]} &= \gamma G_{[3]} + R_2 = R_2 + \gamma R_3 + \gamma^2 R_4 + \gamma^3 R_T, & V[S_N] += G_{[4]}
|
||||
\\\\
|
||||
G_{[3]} &= \gamma G_{[2]} + R_3 = R_3 + \gamma R_4 + \gamma^2 R_T, & V[S_L] += G_{[3]}
|
||||
\\\\
|
||||
G_{[2]} &= \gamma G_{[1]} + R_4 = R_4 + \gamma R_T, & V[S_G] += G_{[2]}
|
||||
\\\\
|
||||
G_{[1]} &= R_T, & V[S_E] += G_{[1]}
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
|
||||
@@ -192,7 +192,7 @@ $n$ 必须大于等于 3，这个游戏才能进行下去。所以当 $n>2$ 时
|
||||
|
||||
### 代码模拟
|
||||
|
||||
代码位置：。。。。。。
|
||||
代码位置：ThreeDoors.py
|
||||
|
||||
如果读者忘记了概率论的基本知识，也不擅长于穷举推导，但是还残留有一点点的编程技巧，那么我们可以用代码模拟上述过程，看看结果如何。这也是程序员学习理论知识的窍门之一。
|
||||
|
||||
|
||||
@@ -41,8 +41,10 @@
|
||||
- 如果是 [...B, C, A...] 的顺序，则 $nB$ 计数，$nA$ 也计数一次。问号表示任意字母。
|
||||
- 如果是 [...B, B, A, A...] 的顺序，则要分别计数两次。
|
||||
|
||||
***代码位置：RentCar_1_FromB.py***
|
||||
|
||||
关键函数代码如下：
|
||||
|
||||
代码片段如下：
|
||||
```Python
|
||||
# 统计rent_from出现的次数，以及经过t天转移到return_to的次数
|
||||
def counter_from_t_to(X, rent_from, return_to, t):
|
||||
@@ -56,7 +58,8 @@ def counter_from_t_to(X, rent_from, return_to, t):
|
||||
return n_from, n_to
|
||||
```
|
||||
|
||||
运行代码 RentCar_1_FromB.py，得到以下结果：
|
||||
运行后得到以下结果：
|
||||
|
||||
```
|
||||
天数 = 2
|
||||
B->A : 2979,238
|
||||
@@ -78,7 +81,9 @@ B->D : 2888,680
|
||||
其实这是一个标准的**转移概率**的问题。
|
||||
|
||||
### 2 转移概率
|
||||
|
||||
|
||||
***代码位置：RentCar_2_OneByOne.py***
|
||||
|
||||
转移概率是马尔可夫链中的重要概念，我们后面再讲马尔科夫链，本节先把转移概率的问题搞清楚。
|
||||
|
||||
在本问题中，有 4 个门店，所以对于某辆车来说，它第二天早晨出现在哪个门店就有 4 种可能，我们称之为 4 个状态：$[A, B, C, D]$。
|
||||
@@ -215,6 +220,8 @@ OK! 2 天后的问题解决了，那么 5 天后呢？这么计算太麻烦了
|
||||
|
||||
### 4 转移概率矩阵
|
||||
|
||||
***代码位置：RentCar_4_Matrix.py***
|
||||
|
||||
**转移概率矩阵**，又称为**状态分布矩阵**，它用矩阵的形式定义了从任意事件（状态）转换到另一个事件（状态）的概率。
|
||||
|
||||
如何得到这个矩阵呢？
|
||||
@@ -319,6 +326,8 @@ $$
|
||||
|
||||
### 5 迭代计算
|
||||
|
||||
***代码位置：RentCar_5_OneByOne.py***
|
||||
|
||||
下面我们使用矩阵计算来代替复杂的循环逻辑：
|
||||
|
||||
1. 首先定义第 0 天的初始向量：
|
||||
@@ -372,7 +381,7 @@ $$
|
||||
|
||||
我们可以很方便地写出代码，来计算第 5 天的出现概率：
|
||||
|
||||
```Python
|
||||
```python
|
||||
P = np.array([
|
||||
[0.1, 0.3, 0.0, 0.6],
|
||||
[0.8, 0.0, 0.2, 0.0],
|
||||
@@ -409,6 +418,8 @@ $[0.32408,\ 0.26466,\ 0.18038,\ 0.23088]$
|
||||
|
||||
### 6 K步转移概率矩阵
|
||||
|
||||
***代码位置：RentCar_6_KStep.py***
|
||||
|
||||
从上面的例子我们可以看到，只有一步的转移概率是不能满足实际需要的，通常需要迭代计算才能知道 K 步后的情况如何，虽然这已经比没有矩阵时的复杂循环逻辑好了很多。
|
||||
|
||||
比如，5 步后的概率应该是：
|
||||
@@ -485,6 +496,8 @@ $$
|
||||
|
||||
### 7 更多步的情况
|
||||
|
||||
***代码位置：RentCar_7_Convergence.py***
|
||||
|
||||
从 K 步可以向更远的地方思考，如果 K 为无穷大时，这个转移概率矩阵是什么情况呢？
|
||||
|
||||
同样可以用代码做一个试验：
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 19 KiB After Width: | Height: | Size: 18 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 12 KiB After Width: | Height: | Size: 12 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 18 KiB |
@@ -1,38 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
# 从酒馆出发运行向家的相反方向走
|
||||
def RandomWalker(distance=50):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
trajectory = [] # 行走的路径
|
||||
while(position < distance):
|
||||
# 随机选择向前(1)向后(-1)不动(0), 概率是[0.4,0.2,0.4]
|
||||
step = np.random.choice([-1,0,1], p=[0.4,0.2,0.4])
|
||||
position += step # 更新位置
|
||||
counter += 1 # 更新步数
|
||||
trajectory.append(position) # 记录位置
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(position, counter, np.mean(trajectory))
|
||||
|
||||
# 从酒馆出发只能向家的方向走
|
||||
def RandomWalker2(distance=50):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
trajectory = [] # 行走的路径
|
||||
while(position < distance):
|
||||
# 随机选择向前(1)向后(-1)不动(0), 概率是[0.4,0.2,0.4]
|
||||
step = np.random.choice([-1,0,1], p=[0.4,0.2,0.4])
|
||||
position += step # 更新位置
|
||||
if (position < 0):
|
||||
position = 0
|
||||
counter += 1 # 更新步数
|
||||
trajectory.append(position) # 记录位置
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(position, counter, np.mean(trajectory))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
for i in range(10):
|
||||
RandomWalker(50)
|
||||
@@ -1,38 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
# 从酒馆出发运行向家的相反方向走
|
||||
def RandomWalker(distance=50):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
trajectory = [] # 行走的路径
|
||||
while(position < distance):
|
||||
# 随机选择向前(1)向后(-1)不动(0), 概率是[0.4,0.2,0.4]
|
||||
step = np.random.choice([-1,0,1], p=[0.4,0.2,0.4])
|
||||
position += step # 更新位置
|
||||
counter += 1 # 更新步数
|
||||
trajectory.append(position) # 记录位置
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(position, counter, np.mean(trajectory))
|
||||
|
||||
# 从酒馆出发只能向家的方向走
|
||||
def RandomWalker2(distance=50):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
trajectory = [] # 行走的路径
|
||||
while(position < distance):
|
||||
# 随机选择向前(1)向后(-1)不动(0), 概率是[0.4,0.2,0.4]
|
||||
step = np.random.choice([-1,0,1], p=[0.4,0.2,0.4])
|
||||
position += step # 更新位置
|
||||
if (position < 0):
|
||||
position = 0
|
||||
counter += 1 # 更新步数
|
||||
trajectory.append(position) # 记录位置
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(position, counter, np.mean(trajectory))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
for i in range(10):
|
||||
RandomWalker(50)
|
||||
@@ -0,0 +1,20 @@
|
||||
import numpy as np
|
||||
|
||||
# 从酒馆出发,允许向家的相反方向走
|
||||
def RandomWalker(distance=50):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
trajectory = [] # 行走的路径
|
||||
while(position < distance):
|
||||
# 随机选择向前(1)向后(-1)不动(0), 概率是[0.4,0.2,0.4]
|
||||
step = np.random.choice([-1,0,1], p=[0.4,0.2,0.4])
|
||||
position += step # 更新位置
|
||||
counter += 1 # 更新步数
|
||||
trajectory.append(position) # 记录位置
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(str.format("步数 : {0}\t最远 : {1}\t平均步数 : {2}", counter, np.min(trajectory), np.mean(trajectory)))
|
||||
|
||||
if __name__ == "__main__":
|
||||
for i in range(10):
|
||||
RandomWalker(50)
|
||||
+2
-17
@@ -1,22 +1,7 @@
|
||||
import numpy as np
|
||||
|
||||
# 从酒馆出发运行向家的相反方向走
|
||||
def RandomWalker(distance=50):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
trajectory = [] # 行走的路径
|
||||
while(position < distance):
|
||||
# 随机选择向前(1)向后(-1)不动(0), 概率是[0.4,0.2,0.4]
|
||||
step = np.random.choice([-1,0,1], p=[0.4,0.2,0.4])
|
||||
position += step # 更新位置
|
||||
counter += 1 # 更新步数
|
||||
trajectory.append(position) # 记录位置
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(position, counter, np.mean(trajectory))
|
||||
|
||||
# 从酒馆出发只能向家的方向走
|
||||
def RandomWalker2(distance=50):
|
||||
def RandomWalker(distance=50):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
trajectory = [] # 行走的路径
|
||||
@@ -30,7 +15,7 @@ def RandomWalker2(distance=50):
|
||||
trajectory.append(position) # 记录位置
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(position, counter, np.mean(trajectory))
|
||||
print(str.format("步数 : {0}\t平均步数 : {1}", counter, np.mean(trajectory)))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
@@ -0,0 +1,23 @@
|
||||
import numpy as np
|
||||
|
||||
# 从酒馆出发,允许向家的相反方向走,但是遇到终点后再也不能回家
|
||||
def RandomWalker(distance_home=50, distance_end=-200):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
trajectory = [] # 行走的路径
|
||||
while(True):
|
||||
# 随机选择向前(1)向后(-1)不动(0), 概率是[0.4,0.2,0.4]
|
||||
step = np.random.choice([-1,0,1], p=[0.4,0.2,0.4])
|
||||
position += step # 更新位置
|
||||
counter += 1 # 更新步数
|
||||
trajectory.append(position) # 记录位置
|
||||
if (position == distance_home or position == distance_end):
|
||||
break
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(str.format("步数 : {0}\t位置 : {1}\t平均步数 : {2}", counter, position, np.mean(trajectory)))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
for i in range(10):
|
||||
RandomWalker(50)
|
||||
@@ -1,5 +1,5 @@
|
||||
|
||||
## 醉汉回家问题 - 马尔可夫理论
|
||||
## 醉汉回家问题 - 马尔可夫过程
|
||||
|
||||
### 1 提出问题
|
||||
|
||||
@@ -23,9 +23,10 @@
|
||||
|
||||
不善于理论推导的话，可以发挥计算机的优势，写一段代码来模拟这个醉汉：
|
||||
|
||||
```Python
|
||||
import numpy as np
|
||||
***代码位置：RandomWalker.py***
|
||||
|
||||
```Python
|
||||
# 从酒馆出发,允许向家的相反方向走
|
||||
def RandomWalker(distance=50):
|
||||
position = 0 # 距离酒馆的位置，为50时表示到家
|
||||
counter = 0 # 行走的步数（包括原地不动）
|
||||
@@ -38,29 +39,29 @@ def RandomWalker(distance=50):
|
||||
trajectory.append(position) # 记录位置
|
||||
|
||||
# 输出最后位置和步数，计算位置平均值
|
||||
print(position, counter, np.mean(trajectory))
|
||||
|
||||
print(str.format("步数 : {0}\t最远 : {1}\t平均步数 : {2}", counter, np.min(trajectory), np.mean(trajectory)))
|
||||
|
||||
if __name__ == "__main__":
|
||||
for i in range(10): # 试验10次
|
||||
RandomWalker(50)
|
||||
```
|
||||
运行上述代码，可以得到 10 次试验结果如下：
|
||||
```
|
||||
50 2383 22.091061686949224
|
||||
50 42234 -59.706681820334325
|
||||
50 11292 -23.148069429684732
|
||||
50 230781 -163.41915062331822
|
||||
50 7055 6.20722891566265
|
||||
50 1366 11.948023426061493
|
||||
50 1387 13.202595529920693
|
||||
50 1742 22.126291618828933
|
||||
50 5477 -5.9565455541354755
|
||||
50 14322 -32.92989805893032
|
||||
步数 : 1404 最远 : -11 平均步数 : 13.004985754985755
|
||||
步数 : 5854 最远 : -40 平均步数 : -4.8963102152374445
|
||||
步数 : 8312 最远 : -84 平均步数 : -15.017204042348412
|
||||
步数 : 7645 最远 : -70 平均步数 : -7.099542184434271
|
||||
步数 : 5420 最远 : -36 平均步数 : 8.072878228782288
|
||||
步数 : 1419 最远 : -4 平均步数 : 15.985200845665961
|
||||
步数 : 22170 最远 : -114 平均步数 : -41.33098782138024
|
||||
步数 : 1699 最远 : -13 平均步数 : 9.114184814596822
|
||||
步数 : 128704 最远 : -276 平均步数 : -90.89893088015913
|
||||
步数 : 1034 最远 : -1 平均步数 : 25.192456479690524
|
||||
```
|
||||
数据解读：
|
||||
|
||||
- 第一次试验结果表面，最终移动了 50 个单元位置到家了，走了 2383 步，历史路径的平均位置是 22。
|
||||
- 第二次试验结果表面，最终移动了 50 个单元位置到家了，走了 42234 步，历史路径的平均位置是 -59.7。
|
||||
- 第一次试验结果，走了 1404 步，最远走到了反向 11 步的地方，历史路径的平均位置是 13。
|
||||
- 第二次试验结果，走了 5854 步，最远走到了反向 40 步的地方，历史路径的平均位置是 -4.89。
|
||||
......
|
||||
|
||||
这个结果说明了醉汉是一定可以到家的，但是步数和平均位置的方差很大。
|
||||
@@ -232,10 +233,11 @@ $$
|
||||
- 处于红色方框的位置已经到家了，处于吸收状态（不再移动），所以到家的概率 $P_0=1$
|
||||
- 中间有 3 个位置，$X+1,X,X-1$，从这三个位置到家的概率分别是 $P_{X+1},P_{X},P_{X-1}$
|
||||
|
||||
<center>
|
||||
<img src="./img/RandomWalker-3.png" width="600">
|
||||
|
||||
图 3
|
||||
|
||||
</center>
|
||||
|
||||
由于从 X-1 到 X 位置有0.5的概率，从 X+1 到 X 位置也有0.5的概率，所以
|
||||
$$
|
||||
@@ -284,12 +286,13 @@ $$
|
||||
从式 8 可以看出，$n$ 越大，醉汉回家的概率越大，$n \to \infty$ 时，醉汉回家的概率接近于 1，这是我们的问题的原意。
|
||||
|
||||
|
||||
随机（相当于醉汉）游走问题是数学史上的一个著名问题，1905年，英国统计学家Pearson在《自然》杂志上公开求解随机游走（Random Walk）问题。1921年，匈牙利数学家波利亚（Polya，1887-1985）在研究随机游走问题后，提出了著名的随机游走定理，证明一维或二维随机游走返回原点的概率为100%，从而得出了醉汉最终会返回原点的结论。Polya随机游走定理被《The Math Book》誉为数学史上250个里程碑式的重大发现之一，Polya本人也被人们视为20世纪最具影响力的数学家之一。日本著名数学家角谷静夫通俗形象地将Polya随机游走定理表述为：喝醉的酒鬼总能找到回家的路。因此，随机游走定理也被称为酒鬼回家定理。
|
||||
随机（相当于醉汉）游走问题是数学史上的一个著名问题，1905年，英国统计学家 Pearson 在《自然》杂志上公开求解随机游走（Random Walk）问题。1921年，匈牙利数学家波利亚（Polya，1887-1985）在研究随机游走问题后，提出了著名的随机游走定理，证明一维或二维随机游走返回原点的概率为 100%，从而得出了醉汉最终会返回原点的结论。Polya 随机游走定理被《The Math Book》誉为数学史上 250 个里程碑式的重大发现之一，Polya 本人也被人们视为20世纪最具影响力的数学家之一。日本著名数学家角谷静夫通俗形象地将 Polya 随机游走定理表述为：喝醉的酒鬼总能找到回家的路。因此，随机游走定理也被称为酒鬼回家定理。
|
||||
|
||||
随机游走是概率论与随机过程学科中用于描述随机现象的一种基本随机过程。液体中悬浮微粒的布朗运动、空气中的烟雾扩散、光纤陀螺的随机游走误差等动态随机现象均可用随机游走模型进行描述。
|
||||
|
||||
波利亚令人吃惊地证明了在维数比2更高的情况下，酒鬼回家的概率大大小于1！比如说，在三维网格中随机游走，最终能回到出发点的概率只有 34%。
|
||||
波利亚令人吃惊地证明了三维及以上的情况下，酒鬼回家的概率大大小于1！比如说，在三维网格中随机游走，最终能回到出发点的概率只有 34%。
|
||||
|
||||
表 2 空间维度与返回原点概率的关系
|
||||
|
||||
|空间维度|返回原点概率|
|
||||
|:-:|-:|
|
||||
@@ -304,6 +307,55 @@ $$
|
||||
|
||||
酒鬼不可能在空中游走，鸟儿的活动空间才是3维的，因此，日本数学家角谷静夫（Shizuo Kakutani，1911–2004）将波利亚定理用一句通俗又十分风趣的语言来总结：喝醉的酒鬼总能找到回家的路，喝醉的小鸟则可能永远也回不了家。
|
||||
|
||||
### 6 更多试验
|
||||
|
||||
#### 只允许向家的方向走
|
||||
|
||||
***代码位置：RandomWalker_Forward.py***
|
||||
|
||||
当醉汉向相反方向走时，我们强制他留在原地不动，即有 $0.4+0.2=0.6$ 的概率留在酒馆。如图 4 所示。
|
||||
|
||||
<center>
|
||||
<img src="./img/RandomWalker-4.png">
|
||||
|
||||
图 4 从酒馆出发有 0.6 的概率留在酒馆
|
||||
</center>
|
||||
|
||||
这种限制下，情况要好很多，因为不能反向走，所以醉汉更容易回到家。10 次试验的运行结果如下：
|
||||
|
||||
```
|
||||
步数 : 1232 平均步数 : 23.3125
|
||||
步数 : 4199 平均步数 : 16.6799237913789
|
||||
步数 : 2229 平均步数 : 24.635711081202334
|
||||
步数 : 954 平均步数 : 23.31970649895178
|
||||
步数 : 3284 平均步数 : 14.660170523751523
|
||||
步数 : 1680 平均步数 : 16.552380952380954
|
||||
步数 : 4645 平均步数 : 16.27491926803014
|
||||
步数 : 1308 平均步数 : 16.874617737003057
|
||||
步数 : 414 平均步数 : 26.08937198067633
|
||||
步数 : 2824 平均步数 : 10.885977337110482
|
||||
```
|
||||
|
||||
步数和平均步数比第一个试验要少很多。
|
||||
|
||||
#### 可以反向走，但是有失败的限制
|
||||
|
||||
为了理解上一小节的理论推导，可以把图 3 中的 n 设置为 250，即，反向到达 200 步的地方，就认为是回家失败（醉汉被野兽吃掉了）。
|
||||
|
||||
```
|
||||
步数 : 9653 位置 : 50 平均步数 : -2.3523257018543458
|
||||
步数 : 24152 位置 : 50 平均步数 : -29.24494865849619
|
||||
步数 : 12230 位置 : -200 平均步数 : -83.99648405560099
|
||||
步数 : 56444 位置 : -200 平均步数 : -83.09604209481964
|
||||
步数 : 1740 位置 : 50 平均步数 : 12.133333333333333
|
||||
步数 : 23286 位置 : 50 平均步数 : -60.63591857768616
|
||||
步数 : 15289 位置 : -200 平均步数 : -74.86349663156518
|
||||
步数 : 24192 位置 : 50 平均步数 : -9.22829861111111
|
||||
步数 : 5952 位置 : 50 平均步数 : -18.80342741935484
|
||||
步数 : 2043 位置 : 50 平均步数 : 12.900636319138522
|
||||
```
|
||||
运行 10 次试验后，可以看到其中有 3 次都失败了，即位置达到了 -200 的地方。
|
||||
|
||||
|
||||
### 参考资料
|
||||
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 16 KiB |
@@ -1,27 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
P = np.array(
|
||||
[ #Game Cl1 Cl2 Cl3 Pass Rest End
|
||||
[0.9, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.2],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.6, 0.4, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.2, 0.4, 0.4, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
|
||||
]
|
||||
)
|
||||
|
||||
def Check_Convergence(P):
|
||||
P_curr = P.copy()
|
||||
for i in range(100000):
|
||||
P_next=np.dot(P,P_curr)
|
||||
print("迭代次数 =",i+1)
|
||||
print(np.around(P_next, 2))
|
||||
if np.allclose(P_curr, P_next, rtol=1e-2, atol=1e-4):
|
||||
break
|
||||
P_curr = P_next
|
||||
return P_next
|
||||
|
||||
if __name__=="__main__":
|
||||
Pn = Check_Convergence(P)
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 37 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 41 KiB |
|
Before Width: | Height: | Size: 26 KiB After Width: | Height: | Size: 26 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 22 KiB |
|
Before Width: | Height: | Size: 13 KiB After Width: | Height: | Size: 13 KiB |
|
Before Width: | Height: | Size: 12 KiB After Width: | Height: | Size: 12 KiB |
|
Before Width: | Height: | Size: 14 KiB After Width: | Height: | Size: 14 KiB |
|
Before Width: | Height: | Size: 12 KiB After Width: | Height: | Size: 12 KiB |
@@ -0,0 +1,46 @@
|
||||
import numpy as np
|
||||
|
||||
P = np.array(
|
||||
[
|
||||
[0.0, 0.9, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.2, 0.1, 0.1, 0.1, 0.3, 0.1, 0.0, 0.1, 0.0],
|
||||
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.7, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.2, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.5, 0.0, 0.0],
|
||||
[0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
|
||||
[0.0, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
|
||||
]
|
||||
)
|
||||
|
||||
# |出发|正常行驶|礼让行人|闹市减速|超速行驶|路口闯灯|小区减速|拨打手机|发生事故|安全抵达|结束|
|
||||
R = [0, 0, +1, +1, -3, -6, +1, -3, -1, +5, 0]
|
||||
|
||||
def Matrix(gamma):
|
||||
num_state = P.shape[0]
|
||||
I = np.eye(len(R)) * (1+1e-7)
|
||||
tmp1 = I - gamma * P
|
||||
tmp2 = np.linalg.inv(tmp1)
|
||||
vs = np.dot(tmp2, R)
|
||||
return vs
|
||||
|
||||
def Check_Convergence(P):
|
||||
P_curr = P.copy()
|
||||
for i in range(100000):
|
||||
P_next=np.dot(P,P_curr)
|
||||
print("迭代次数 =",i+1)
|
||||
print(np.around(P_next, 2))
|
||||
if np.allclose(P_curr, P_next, rtol=1e-2, atol=1e-4):
|
||||
break
|
||||
P_curr = P_next
|
||||
return P_next
|
||||
|
||||
if __name__=="__main__":
|
||||
#Pn = Check_Convergence(P)
|
||||
#print(Matrix(0.9))
|
||||
#print(Matrix(0))
|
||||
print(np.around(Matrix(1),1))
|
||||
|
||||
+138
-104
@@ -1,9 +1,9 @@
|
||||
|
||||
## 学生学习问题
|
||||
## 马尔可夫奖励过程
|
||||
|
||||
### 1 提出问题
|
||||
|
||||
前面学习了马尔可夫链。如何用这种形式来表现、研究在校学生在一门课程中的学习、考试等一些列状态呢？
|
||||
前面学习了马尔可夫链，它可以帮助分析复杂的多状态转移问题。如何用这种形式来表现和研究安全驾驶中的一系列问题呢？
|
||||
|
||||
首先我们提出一个概念：有终止状态的马尔可夫链。
|
||||
|
||||
@@ -17,87 +17,96 @@
|
||||
|
||||
当然，真实世界的马尔科夫链的状态空间很可能是成千上万个，转移路径也更复杂。
|
||||
|
||||
|
||||
|
||||
|
||||
### 2 建立模型
|
||||
|
||||
以在校大学生为例，假设一门课只需要上三次课就可以结束，然后就可以通过考试而结课拿学分。当然，这其中也不是那么顺利的，学生可能会遇到各种挑战：
|
||||
以一个司机驾车上路为例，我们可根据日常驾驶经验以及在路上遇到的各种路况来建立一个模型。图 1 是一个有关安全驾驶的一系列状态的马尔可夫链，也可以叫做状态转移图。
|
||||
|
||||
- 上课不专心听讲而去打手机游戏；
|
||||
- 学到一半的时候觉得这门课索然无味，中途退课；
|
||||
- 觉得离考试还远，不着急复习巩固知识，而是去休息；
|
||||
......
|
||||
<center>
|
||||
<img src="./img/Drive-1.png" width="600">
|
||||
|
||||
图 1 是一个有关学生的学习、考试等一些列状态的马尔可夫链，也可以叫做状态转移图。
|
||||
图 1 安全驾驶问题的状态转移概率图
|
||||
</center>
|
||||
|
||||
<img src="./img/Student-1.png" width="500">
|
||||
状态说明：
|
||||
|
||||
图 1 学习问题的状态转移概率图
|
||||
|
||||
- Class 1,2,3
|
||||
上课/学习/复习状态，假设一门课是需要三次课的学习就可以结束。
|
||||
- 在上课 1 中，有 0.5 的概率跑到 Game 状态，在课堂上用手机偷偷摸摸打游戏，另外 0.5 的概率到上课 2 状态。
|
||||
- 在上课 2 中，有 0.2 的概率直接退课，另外 0.8 的概率到上课 3 状态。
|
||||
- 在上课 3 中，有 0.6 的概率去考试，另外 0.4 的概率以为考试还远，不着急准备呢，就去休息了。
|
||||
- Pass
|
||||
考试通过状态。考试结束后，以 100% 概率到结课状态。
|
||||
- Rest
|
||||
休息状态。休息结束后，发现学的内容都忘得差不多了，遂分别以不同的概率回到三次课的上课/学习/复习状态。
|
||||
在这里我们不区分上课和复习，如果是第一次到达 C1 状态，就认为是上课，否则就认为是复习。
|
||||
- Game
|
||||
娱乐/打游戏状态。在游戏状态中很大可能不能自拔，以 0.9 概率继续打游戏，只有 0.1 的概率幡然悔悟回到学习状态。
|
||||
- End
|
||||
结课/退课状态。进入此状态后将不再进行转移，或者是说以 100% 的概率转移到自己，叫做结束状态或者吸收状态。
|
||||
|
||||
表 1 中列出了状态转移矩阵，与租车问题中的矩阵形式相同。
|
||||
|
||||
表 1 状态转移矩阵
|
||||
|
||||
|P: 从$\rightarrow$到|Game|Class1|Class2|Class3|Pass|Rest|End|
|
||||
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|
||||
|**Game**|0.9|0.1||||||
|
||||
|**Class1**|0.5||0.5|||||
|
||||
|**Class2**||||0.8|||0.2|
|
||||
|**Class3**|||||0.6|0.4||
|
||||
|**Pass**|||||||1.0|
|
||||
|**Rest**||0.2|0.4|0.4||||
|
||||
|**End**|||||||1.0|
|
||||
|
||||
有的读者可能有个疑问：打游戏上瘾，从 Game 到 Game 有 0.9 的高概率，那么当该学生打游戏 2 小时后良心发现，转到学习状态的概率会不会大于 0.1 呢？
|
||||
|
||||
这是一个简单的平稳环境的马尔科夫链，如果考虑更复杂的情况，可以在 Game 状态下增加一个计数器：
|
||||
- 如果打游戏超过 1 小时了，则有 0.5 的概率回到学习状态；
|
||||
- 如果没超过 1 小时，则有 0.1 的概率回到学习状态。
|
||||
|
||||
对于这种有终止状态的状态转移概率图（矩阵），是没有平稳、收敛的概念的，感兴趣的读者可以运行 Studnet.py 来验证结果。
|
||||
- 出发
|
||||
- 0.9 的概率心情较好，进入正常行驶状态；
|
||||
- 0.1 的概率可能有急事而超速行驶。
|
||||
- 正常行驶
|
||||
正常行驶是标准状态，但是由于各种路况，可能会转移到其它几个状态：
|
||||
- 0.1 的概率由于路况好或与他人斗气飙车而超速行驶；
|
||||
- 0.1 的概率开车时心不在焉，在路口闯灯；
|
||||
- 0.2 的概率在斑马线前礼让行人；
|
||||
- 0.1 的概率开到闹市时减速行驶；
|
||||
- 0.1 的概率遇到急事而拨打手机；
|
||||
- 0.4 的概率进入目的地区域后减速行驶。
|
||||
- 拨打手机
|
||||
指的是开车过程中拨打手机，属于危险行为。
|
||||
- 0.6 的概率结束通话返回正常行驶；
|
||||
- 0.4 的概率发生事故。
|
||||
- 礼让行人
|
||||
礼让后以 1.0 的概率回到正常行驶状态。
|
||||
- 小区减速
|
||||
减速后以 1.0 的概率安全抵达终点。
|
||||
- 安全抵达
|
||||
安全抵达后进入结束状态。
|
||||
- 闹市减速
|
||||
在闹市低速行驶时：
|
||||
- 0.3 的概率遇到行人较多，礼让；
|
||||
- 0.7 的概率回到正常行驶状态。
|
||||
- 超速行驶
|
||||
属于危险行为
|
||||
- 0.2 的概率回到正常行驶；
|
||||
- 0.3 的概率闯红灯；
|
||||
- 0.5 的概率发生事故。
|
||||
- 路口闯灯
|
||||
- 0.9 的概率会出事故；
|
||||
- 0.1 的概率遇到警察心情好，侥幸回到正常行驶。
|
||||
- 发生事故
|
||||
1.0 的概率结束，不能再达到目的地。
|
||||
- 结束
|
||||
进入此状态后将不再进行转移，或者是说以 100% 的概率转移到自己，叫做结束状态或者吸收状态。
|
||||
|
||||
### 3 分幕和采样
|
||||
|
||||
由于终止状态的存在（图 1 中的 End 状态），可以引入一个新的强化学习中的重要概念：**分幕**（Episode）。比如：
|
||||
|
||||
- 在醉汉回家问题中，醉汉到家了，整个过程叫做一幕；第二天该醉汉又喝醉了，再次到家后，又叫做一幕。
|
||||
- 在学生学习问题中，学生经过一系列折腾，最终到达 End 状态，叫做一幕。
|
||||
- 在安全驾驶问题中，司机经过一系列复杂的路况，最终到达目的地或者出事故，叫做一幕。
|
||||
- 一盘棋中，最终双方分出输赢或者打平，叫做一幕。
|
||||
- 打桥牌中，每个人出完手中的最后一张牌，桌面上有 13 墩牌，叫做一幕。
|
||||
- 一个扫地机器人，完成当天任务回到充电状态，叫做一幕；如果在打扫过程中电量过低，不得不暂时回到充电状态，不叫作完成一幕。
|
||||
|
||||
在图 1 中，可以根据不同的学生在本门课中的经历，获得不同的到达终点路径，其中，用 C1,C2,C3 表示 Class 1,Class 2,Class 3。比如：
|
||||
- C1 - C2 - C3 - Pass - End
|
||||
- C1 - Game - Game - Game - C1 - C2 - End
|
||||
- C1 - C2 - C3 - Rest - C2 - C3 - Pass - End
|
||||
- C1 - C2 - C3 - Rest - C1 - Game - Game - C1 - C2 - End
|
||||
在图 1 中，可以根据不同的司机上路的经历，获得不同的到达终点路径，比如：
|
||||
- S:出发 - N:正常行驶 - L:小区减速 - G:安全抵达 - E:结束
|
||||
- S:出发 - N:正常行驶 - P:礼让行人 - N:正常行驶 - L:小区减速 - G:安全抵达 - E:结束
|
||||
- S:出发 - N:正常行驶 - R:路口闯灯 - C:发生事故 - E:结束
|
||||
- S:出发 - X:超速行驶 - R:路口闯灯 - N:正常行驶 - R:路口闯灯 - C:发生事故 - E:结束
|
||||
......
|
||||
|
||||
上述的过程叫做**采样**（Sample）。对于大多数学生来说，经历的都是第一个路径，而有极少数学生走的第四个路径。
|
||||
上述的过程叫做**采样**（Sample）。对于大多司机来说，经历的都是第一个路径，而有极少司机走的第四个路径（估计是酒后驾车）。
|
||||
|
||||
根据转移概率，经过采样后到达终点，得到一幕状态序列。没有到达终点的状态序列不叫作完整的状态序列，也不叫做一幕，但是这些序列片段仍然有研究价值。
|
||||
|
||||
另外，读者可能注意到在上述采样的例子中，都是从 C1 开始的，似乎有一个隐含的开始状态。在某些问题中，我们指定一些状态为开始状态，原因是：
|
||||
另外，读者可能注意到在上述采样的例子中，都是从 S 开始的，似乎有一个隐含的开始状态。在某些问题中，我们指定一些状态为开始状态，原因是：
|
||||
|
||||
- 一是为了符合实际问题的逻辑需要，比如学生学习问题，一般都是从第一节课开始的。
|
||||
- 一是为了符合实际问题的逻辑需要，比如安全驾驶问题，一般都是从出门开始的。
|
||||
- 二是为了不过分简化问题的难度。比如一个迷宫，如果不指定开始状态，而是从快到出口的一个位置作为起始位置，将会大大降低问题的难度。
|
||||
|
||||
### 4 奖励和收益
|
||||
|
||||
从图 1 中，读者可以体会到，在 C3 状态时，会比在 C1 状态有更多的可能性到达考试通过并结课的状态的（不考虑退课的情况），那么 C3 状态应该比 C1 状态要“好”。如何定义这个“好”呢？这就要引入强化学习的两个重要概念：**奖励**（Reward），以及**回报**（Return 或 Gain）。
|
||||
通过学习交通法规以及上路实践，读者会知道：
|
||||
- 如果超速行驶，会面临至少 3 分的扣分；
|
||||
- 如果闯红灯，扣 6 分；
|
||||
- 如果驾驶时拨打手机，扣 3 分
|
||||
......
|
||||
|
||||
在自动驾驶中，如何让智能体也“懂得”这些交通规则呢？
|
||||
|
||||
一个办法是制定一些规则，以代码的形式硬编到逻辑中；另外一个办法是通过学习，知道什么可以做，什么不可以做。这就要引入强化学习的两个重要概念：**奖励**（Reward），以及**回报**（Return）。
|
||||
|
||||
#### 奖励
|
||||
|
||||
@@ -107,25 +116,39 @@
|
||||
|
||||
奖励有两种定义方式，如图 2 所示。
|
||||
|
||||
<center>
|
||||
<img src="./img/Reward-1.png" width="600">
|
||||
|
||||
图 2 奖励的两种定义方式
|
||||
</center>
|
||||
|
||||
- 给状态定义奖励
|
||||
图 2 左图中，从状态 $S_a$ 或 $S_b$ 到达状态 $S_c$ 时，获得同样的奖励 $R$，没有差别。
|
||||
比如，走一个迷宫，别人找到了最佳路径用了 1 分钟完成任务，而你走了弯路，3 分钟才完成任务，但是裁判并不会因为你花费了更多的力气而给你更多的奖励，**注重结果**。
|
||||
- 给状态定义奖励 —— **注重结果**
|
||||
图 2 左图中：
|
||||
- 从状态 $S_a$ 到达 $S_c$ 后，获得奖励 $R_1$。
|
||||
- 从状态 $S_a$ 到达 $S_d$ 后，获得奖励 $R_2$。
|
||||
- 从状态 $S_b$ 到达 $S_d$ 后，获得奖励 $R_2$。
|
||||
|
||||
可以看到，只要到达 $S_d$，就可以获得 $R_2$，**奖励是给与目标状态的，与源状态无关**。比如，走一个迷宫，别人找到了最佳路径用了 1 分钟完成任务，而你走了弯路，3 分钟才完成任务，但是裁判并不会因为你花费了更多的力气而给你更多的奖励，**注重结果**。
|
||||
|
||||
- 给过程定义奖励
|
||||
图 2 右图中，从状态 $S_a$ 到达状态 $S_c$ 时，获得奖励 $R_1$；从状态 $S_b$ 到达状态 $S_c$ 时，获得奖励 $R_2$。
|
||||
比如，同样是考上清华大学，一个城市里的学生和一个大山里的学生所经历的历程不一样，受教育的环境也不同，应该获得不同的奖励值，**注重过程**。
|
||||
- 给过程定义奖励 —— **注重过程**
|
||||
图 2 右图中：
|
||||
- 从状态 $S_a$ 到达 $S_c$ 时，获得奖励 $R_1$。
|
||||
- 从状态 $S_a$ 到达 $S_c$ 时，获得奖励 $R_2$。
|
||||
- 从状态 $S_b$ 到达 $S_c$ 时，获得奖励 $R_3$。
|
||||
|
||||
注意，这个**奖励与源状态和目标状态的组合有关**。比如，同样是考上清华大学，一个城市里的学生和一个大山里的学生所经历的历程不一样，受教育的环境也不同，应该获得不同的奖励值，**注重过程**。
|
||||
|
||||
更正式的定义是图 2 右图中的方法，即**注重过程**，用数学语言表述就是 $S \times S' \to R$，意为奖励 $R$ 由 $S,S'$ 共同决定，不同的过程会有不同的奖励。但是有时候为了研究问题方便，我们会使用图 2 左图中的定义。这两种方式很容易区分，读者只需要看 $R$ 在图中标注的位置就可以了。
|
||||
更正式的定义是图 2 右图中的方法，即**注重过程**，用数学语言表述就是 $S \times S' \to R$，意为奖励 $R$ 由 $S,S'$ 共同决定，不同的过程会有不同的奖励。但是在目前阶段为了研究问题方便，我们会使用图 2 左图中的定义，即 $S' \to R$。这两种方式很容易区分，读者只需要看 $R$ 在图中标注的位置就可以了。
|
||||
|
||||
有没有一种方式可以统一这两种定义方式呢？在后面学习贝尔曼方程时，再来具体解释。在那之前，我们一直会用**注重结果**的方式，因为它对初学者来说比较容易理解，计算也简便。
|
||||
|
||||
一个完整的奖励过程如图 3 所示。就是在状态转移图中，给每个状态都定义一个奖励值，当到达这个状态时，强化学习过程就会获得相应的奖励值，使得整个过程向着获得最大收益的方向优化和运行。
|
||||
|
||||
<img src="./img/Reward-2.png" width="600">
|
||||
<center>
|
||||
<img src="./img/Reward-2.png">
|
||||
|
||||
图 3 状态与奖励
|
||||
</center>
|
||||
|
||||
很明显，这是用**注重结果**的方式来定义奖励。如果是**注重过程**的话，图 3 的 $S_1$ 状态应该没有奖励值，因为看上去它似乎是起始状态，没有任何**过程**可以定义它的奖励。
|
||||
|
||||
@@ -137,9 +160,9 @@
|
||||
......
|
||||
- 到达终点 $T$ 时，会得到$R_{T}$的奖励。
|
||||
|
||||
这些状态的下标值，只表示前后发生的顺序，即**时刻**，而与状态的序号无关。比如一个状态集中有 4 个状态 $[S_a,S_b,S_c,S_T]$，它们发生的顺序有可能是 $S_a,S_b,S_a,S_c,S_T$，那么有：$S_1=S_a,\ S_2=S_b,\ S_3=S_a,\ S_4=S_c,\ S_5=S_T$。
|
||||
这些状态的下标值，只表示前后发生的顺序，即**时刻**，而与状态的序号无关。比如一个状态集中有 4 个状态 $[S_a,S_b,S_c,S_T]$，马尔可夫链的顺序有可能是 $S_a,S_b,S_a,S_c,S_T$，那么有：$S_1=S_a,\ S_2=S_b,\ S_3=S_a,\ S_4=S_c,\ S_5=S_T$。
|
||||
|
||||
*注：在本书中使用这种定义方式：整个序列是 $S_1,R_2,S_2,R_3,\cdots,S_t,R_{t+1},\cdots$ 的过程。而在有些资料中，用这种定义方式：$S_1,R_1,S_2,R_2,\cdots,S_t,R_{t},\cdots$，需要读者事先注意。*
|
||||
*注：在本书中使用这种定义方式：整个序列的序号是 $S_1,R_2,S_2,R_3,\cdots,S_t,R_{t+1},\cdots$ 的过程。而在有些资料中，用这种定义方式：$S_1,R_1,S_2,R_2,\cdots,S_t,R_{t},\cdots$，需要读者事先注意。*
|
||||
|
||||
#### 回报
|
||||
|
||||
@@ -156,9 +179,11 @@ $$
|
||||
|
||||
智能体的目标是最大化其收到的总收益，即回报。这意味着需要最大化的**不是当前收益，而是长期的累积收益**。我们所有的“目标”或“目的”都可以归结为：最大化智能体接收到的标量信号累积和的概率期望值。
|
||||
|
||||
<img src="./img/Gain.png" width="600">
|
||||
<center>
|
||||
<img src="./img/Gain.png">
|
||||
|
||||
图 4 奖励与回报
|
||||
</center>
|
||||
|
||||
图 4 中展示了 $G_1$ 和 $G_t$ 的计算方式，表示了 $S_1$ 和 $S_t$ 的回报，同时也告诉读者，在一个完整的序列中，我们可以从任意时刻 $t$ 开始计算 $S_t$ 状态的回报，而忽略前面的数据。
|
||||
|
||||
@@ -170,13 +195,31 @@ $$
|
||||
|
||||
### 5 马尔可夫奖励过程（Markov Reward Process）
|
||||
|
||||
用一个文字公式来表示 MRP（Markov Reward Process，马尔可夫奖励过程）：
|
||||
|
||||
$$
|
||||
马尔可夫奖励过程 = 马尔可夫链 + 奖励函数
|
||||
$$
|
||||
|
||||
根据上面学习的知识，再结合上面学生学习的状态转移图 1，我们给每个状态定义一个奖励，如图 5 所示。
|
||||
|
||||
<img src="./img/Student-2.png" width="500">
|
||||
<center>
|
||||
<img src="./img/Drive-2.png" width="600">
|
||||
|
||||
图 5 学生学习问题马尔可夫奖励过程
|
||||
图 5 安全驾驶问题的马尔可夫奖励过程
|
||||
</center>
|
||||
|
||||
显然，我们使用了图 2 中的第一种方式（**注重结果**）来定义奖励，举例来说，无论状态 Class 2 是通过什么路径到达的（可以通过 Class 1 到达，也可以通过 Rest 到达），都可以得到 -2 的奖励。
|
||||
显然，我们使用了图 2 中的第一种方式（**注重结果**）来定义奖励，举例来说，无论状态“发生事故”是通过什么路径到达的（可以通过“拨打手机”到达，也可以通过“路口闯灯”到达），都可以得到 -12 的奖励（实际上是惩罚）。
|
||||
|
||||
奖励函数（值）的设计一般是人工设定的，是通过分析目标问题的实际科学意义或者人文意义来决定的。比如，在图 5 中，通过交通规则的学习，给出制定奖励的过程如下：
|
||||
|
||||
1. 按交规，超速行驶扣 3 分，开车打电话扣 3 分，闯红灯扣 6 分。
|
||||
2. 发生事故扣 1 分。有的读者会有疑问：为什么出了事故只扣 1 分？因为在交规中，除了事故后，不会因为出事故本身而扣分，而是分析出事故的原因，对原因扣分。所以，这里只是象征性地扣 1 分。
|
||||
3. 礼让行人在交规上不加分，但是在强化学习系统中可以加 1 分，以鼓励自动驾驶的智能体强化此状态，保证安全。
|
||||
4. 正常行驶是一个常见状态，得 0 分；
|
||||
5. 闹市减速、小区内减速，和礼让行人一样，都给 1 分奖励。
|
||||
6. 安全抵达给 5 分奖励。
|
||||
7. 出发和结束都是 0 分。
|
||||
|
||||
根据状态转移可以得到一些完整的分幕采样，从而可以计算出每个采样的回报值，列在表 1 中。
|
||||
|
||||
@@ -184,15 +227,12 @@ $$
|
||||
|
||||
||分幕采样序列|回报值计算|
|
||||
|-|-|-|
|
||||
|1|C1-C2-C3-Pass-End|$G_{C1}=-2-2-3+10+0=3$|
|
||||
|2|C1-Game-Game-Game-C1-C2-End|$G_{C1}=-2-1-1-1-2-2+0=-9$|
|
||||
|3|C2-End|$G_{C2}=-2+0=-2$|
|
||||
|4|C1-C2-C3-Rest-C2-C3-Pass-End|$G_{C1}=-2-2-3+1-2-2+10+0=0$|
|
||||
|5|C1-C2-C3-Rest-C1-Game-Game-C1-C2-End|$G_{C1}=-2-2-2+1-2-1-1-2-2+0=-13$|
|
||||
|6|C3-Rest-C1-Game-Game-C1-C2-End|$G_{C3}=-2+1-2-1-1-2-2+0=-9$|
|
||||
|7|Game-Game-C1-C2-End|$G_{Game}=-1-1-2-2+0=-6$|
|
||||
|1|Start - N - L - G - End|$G_{S}=0+0+1+5+0=6$|
|
||||
|2|Start - N - P - N - L - G - End|$G_{S}=0+0+1+0+1+5+0=7$|
|
||||
|3|Start - N - R - C - End|$G_{S}=0+0-6-1+0=-7$|
|
||||
|4|Start - X - R - N - R - C - End|$G_{S}=0-3-6+0-6-1+0=-16$|
|
||||
|
||||
读者可能会产生怀疑：为什么表 1 中的第 1,2,4,5 行都是同样计算 $G_{C1}$ 的回报值，但是有不同的累积结果？这是因为采样不同，路径不同，造成的回报值不同，这种情况是正常的。
|
||||
读者可能会产生怀疑：为什么表 1 中都是同样计算 $G_{S}$ 的回报值，但是有不同的结果？这是因为采样不同，路径不同，造成的回报值不同，这种情况是正常的。
|
||||
|
||||
到目前可以总结出，马尔可夫奖励过程是一个元组的数据序列：$<S,P,R>$，分别表示状态 $S$、转移概率 $P$、奖励 $R$。
|
||||
|
||||
@@ -200,9 +240,7 @@ $$
|
||||
|
||||
我们学习了什么是奖励，但是奖励函数是什么？
|
||||
|
||||
函数这个词可以很宽泛，不一定非得用数学表达式才能表达出来，也可以不是连续的。比如，我们可以给状态 $[S_1, S_2, S3]$ 定义奖励为 $[0, -1, 5]$，这也可以称为“函数”，所以它只是一种“定义”方法。
|
||||
|
||||
比如图 5，可以给状态 $[Class1,Class2,Class3,Pass,Rest,Game,End]$ 定义一个向量作为奖励函数：$[-2,-2,-2,\ 10,\ 1,-1,\ 0]$。
|
||||
函数这个词可以很宽泛，不一定非得用数学表达式才能表达出来，也可以不是连续的。比如图 5，，我们可以给状态 $[S_1, S_2, S3]$ 定义奖励为 $[0, -1, 5]$，这也可以称为“函数”，所以它只是一种“定义”方法。
|
||||
|
||||
而在一些复杂的问题中，确实需要奖励函数，而非简单的奖励值，最常见的有：
|
||||
- 稀疏奖励
|
||||
@@ -221,31 +259,37 @@ $$
|
||||
|
||||
举一个简单的例子来初步理解奖励函数设计的概念。
|
||||
|
||||
<center>
|
||||
<img src="./img/Reward-3.png" width="600">
|
||||
|
||||
图 5 奖励函数的设计
|
||||
图 6 奖励函数的设计
|
||||
</center>
|
||||
|
||||
如图 5 中所示，状态 $S_1,S_2,S_3$ 在相互转移时都没有奖励，或者奖励同为 0，只有到状态 $T$ 时才有 +100 的奖励。为了加快学习进度，把奖励函数做如图 6 的修改。
|
||||
如图 6 中所示，状态 $S_1,S_2,S_3$ 在相互转移时都没有奖励，或者奖励同为 0，只有到状态 $T$ 时才有 +100 的奖励。为了加快学习进度，把奖励函数做如图 7 的修改。
|
||||
|
||||
<center>
|
||||
<img src="./img/Reward-4.png" width="600">
|
||||
|
||||
图 6 奖励函数的修改
|
||||
图 7 奖励函数的修改
|
||||
</center>
|
||||
|
||||
图 6 的学习难度降低了，速度也加快了，但是学习过程很可能陷入 $S_1,S_2,S_1,S_2,\cdots$的循环中（如红色箭头所示），循环了 100+ 次以后，可以得到比 +100 更高的回报，根本不用到达目标状态 $S_T$，这就违背了我们的初衷。
|
||||
图 7 的学习难度降低了，速度也加快了，但是学习过程很可能陷入 $S_1,S_2,S_1,S_2,\cdots$的循环中（如红色箭头所示），循环了 100+ 次以后，可以得到比 +100 更高的回报，根本不用到达目标状态 $S_T$，这就违背了我们的初衷。
|
||||
|
||||
特别地，奖励信号并不是传授智能体如何实现目标的先验知识。例如：
|
||||
- 国际象棋智能体只有当最终获胜时才能获得奖励，而并非达到某个子目标，比如吃掉对方的子或者控制中心区域。如果实现这些子目标也能得到奖励，那么智能体可能会找到某种即使绕开最终目的也能实现这些子目标的方式。例如它可能会找到一种以输掉比赛为代价的方式来吃对方的子。奖励信号只能用来传达什么是你想要实现的目标，而不是如何实现这个目标。
|
||||
- 在狼吃羊的强化学习训练中，如果给狼的每一步移动的奖励设置为 -1，吃到羊的奖励为 +100，那么狼需要在 100 步之内吃到羊才会有正的回报，而在大多数情况下是负值。所以狼可能会选择一开始就一头撞死在障碍物上，以获得 0 分的回报。
|
||||
|
||||
如何避免上述情况呢？我们可以做如图 7 的修改，使得智能体在两个状态间循环时只能得到回报为 0 的过程。
|
||||
如何避免上述情况呢？我们可以做如图 8 的修改，使得智能体在两个状态间循环时只能得到回报为 0 的过程。
|
||||
|
||||
<center>
|
||||
<img src="./img/Reward-5.png" width="600">
|
||||
|
||||
图 7 最终的奖励函数
|
||||
图 8 最终的奖励函数
|
||||
</center>
|
||||
|
||||
### 7 折扣因子 $\gamma$
|
||||
|
||||
如果把**回报**定义为**奖励**的简单相加的话，整个学习框架就会失去一些“灵动”，没有可以调节收益信号大小的“开关”，甚至带来如图 6 所示的灾难。
|
||||
如果把**回报**定义为**奖励**的简单相加的话，整个学习框架就会失去一些“灵动”，没有可以调节收益信号大小的“开关”，甚至带来如图 7 所示的灾难。
|
||||
|
||||
如何避免这个问题呢？
|
||||
|
||||
@@ -282,23 +326,13 @@ $$
|
||||
|
||||
||分幕采样序列|回报值计算（$\gamma=0.9$）|
|
||||
|-|-|-|
|
||||
|1|C1-C2-C3-Pass-End|$G_{C1}=(-2)+0.9*(-2)+0.9^2*(-2)+0.9^3*10+0.9^4*0=1.87$|
|
||||
|2|C1-Game-Game-Game-C1-C2-End|$G_{C1}=(-2)+0.9*(-1)+0.9^2*(-1)+0.9^3*(-1)+0.9^4*(-2)+0.9^5*(-2)+0.9^6*0=-6.93$|
|
||||
|3|C2-End|$G_{C2}=-2+0.9*0=-2$|
|
||||
|...|...|...|
|
||||
|7|Game-Game-C1-C2-End|$G_{Game}=(-1)+0.9*(-1)+0.9^2*(-2)+0.9^3*(-2)+0.9^4*0=-4.98$|
|
||||
|1|Start - N - L - G - End|$G_{S}=0+0.9*0+0.9^2*1+0.9^3*10+0.9^4*0=8.1$|
|
||||
|2|Start - N - P - N - L - G - End|$G_{S}=0+0.9*0+0.9^2*1+0.9^3*0+0.9^4*1+0.9^5*10+0.9^6*0=7.371$|
|
||||
|3|Start - N - R - C - End|$G_{S}=0+0.9*0-0.9^2*6-0.9^3*12+0.9^4*0=-13.608$|
|
||||
|4|Start - X - R - N - R - C - End|$G_{S}=0-0.9*3-0.9^2*6+0.9^3*0-0.9^4*6-0.9^5*12+0.9^6*0=-18.58$|
|
||||
|
||||
表 3 折扣为 0 的分幕采样和回报计算
|
||||
|
||||
||分幕采样序列|回报值计算（$\gamma=0$）|
|
||||
|-|-|-|
|
||||
|1|C1-C2-C3-Pass-End|$G_{C1}=-2$|
|
||||
|2|C1-Game-Game-Game-C1-C2-End|$G_{C1}=-2$|
|
||||
|3|C2-End|$G_{C2}=-2$|
|
||||
|...|...|...|
|
||||
|7|Game-Game-C1-C2-End|$G_{Game}=-1$|
|
||||
|
||||
表 3 中，如果折扣为 0，则回报值 $G$ 就等于当前状态的奖励值 $R$，即$G_t = R_{t+1}$。
|
||||
如果折扣为 0，则回报值 $G$ 就等于当前状态的奖励值 $R$，即 $G_t = R_{t+1}$。
|
||||
|
||||
在式 2 中，如果 $T$ 很大的话，似乎 $G$ 就会变得很大，从而无法计算。但由于限定 $\gamma \le 1$，所以回报值还是不会大得离谱的。特别低，如果奖励值为常数 +1，则回报是：
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 34 KiB |
@@ -1,345 +0,0 @@
|
||||
## 学生问题
|
||||
|
||||
### 1 提出问题
|
||||
|
||||
在学生学习的模型中有很多状态，如何确定某个状态比另一个状态好呢？或者说如何比较两个状态的好坏呢？因为从直觉上讲，学生在 C1,C2,C3 的状态明显要比 Game 状态好，但是如何能用数值的大小来体现这种好坏关系呢？
|
||||
|
||||
<img src="./img/Student-2.png" width="500">
|
||||
|
||||
图 1 学生学习模型
|
||||
|
||||
在上一节中，已经有了分幕、奖励、回报的概念，这一节中，将会利用这些基础概念来定义每个**状态价值函数**，从而可以比较状态之间的好坏。
|
||||
|
||||
需要再次说明的是，在图 1 中，我们使用了**注重结果**的奖励定义方式，直接给每个状态赋值一个奖励，意味只要达到这个状态，就可以立刻获得标注出的奖励值，而不管是从哪条路径达到的。
|
||||
|
||||
### 2 状态价值函数（State Value Function）
|
||||
|
||||
在上一节通过对 $G$ 的计算，以及对状态图的分析理解，我们似乎已经得到了一些启示：距离终点越近的状态，越接近于成功，它的状态价值就越高，似乎用 $G$ 值就可以表示该时刻的状态价值。
|
||||
|
||||
但是，会有一个麻烦出现：每个学生走的路径都不完全一样，状态图虽然是有向的，但是由于环状转移的原因，一个状态在一个序列中可能会出现多次。比如 C1 状态，有三种情况可以到达：
|
||||
|
||||
1. 最开始上第一次课时；
|
||||
2. 从游戏状态返回时；
|
||||
3. 从休息状态返回时。
|
||||
|
||||
而以 C1 开始的路径又可以有很多种，如表 1 所示。
|
||||
|
||||
表 1 以 C1 开始的 $G$ 值的计算
|
||||
|
||||
||分幕采样序列|回报值计算（$\gamma=1$）|
|
||||
|-|-|-|
|
||||
|1|C1-C2-C3-Pass-End|$G_{C1}=-2-2-3+10+0=3$|
|
||||
|2|C1-Game-Game-C1-C2-End|$G_{C1}=-2-1-1-1-2-2+0=-9$|
|
||||
|3|C1-C2-C3-Rest-C2-C3-Pass-End|$G_{C1}=-2-2-3+1-2-2+10+0=0$|
|
||||
|4|C1-C2-C3-Rest-C1-Game-Game-C1-C2-End|$G_{C1}=-2-2-2+1-2-1-1-2-2+0=-13$|
|
||||
|
||||
这样的话，一个状态 C1 就可能有 4 个 $G$ 值，我们用哪个当作其价值函数呢？另外，仔细观察 $G$ 的表达式，它只与时刻及奖励有关，没有体现出状态来。
|
||||
|
||||
考虑到以上两点，定义状态价值函数如下：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
V_t(s) &= \mathbb E [G_t | S_t = s]
|
||||
\\\\
|
||||
&=\mathbb E [ R_{t+1}+\gamma R_{t+2}+\gamma^2 R_{t+3}+ \gamma^3 R_{t+4}+ \cdots]
|
||||
\end{aligned}
|
||||
\tag{1}
|
||||
$$
|
||||
|
||||
### 3 数学期望
|
||||
|
||||
简单地回忆一下**数学期望**的概念。
|
||||
|
||||
一个正常的六面的骰子，投出去后可以得到 [1,2,3,4,5,6] 六种结果，而且概率相等，那么这个骰子的期望值是 $(1+2+3+4+5+6)/6=3.5$。哈哈，读者可能会发现 3.5 这个数子，骰子无法投出来，所以它只是一种定义。
|
||||
|
||||
但是，一个不正常的骰子，比如 [4,5,6] 出现的概率 $p$ 都是 $\frac{1}{5}$，而 [1,2,3] 出现的概率 $p$ 都是 $\frac{2}{15}$，那么它的数学期望是：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
\mathbb E[骰子]&=\sum_{i=1}^6 p_i V_i
|
||||
\\\\
|
||||
&= \frac{2}{15} \times 1+\frac{2}{15} \times 2+\frac{2}{15} \times 3+\frac{1}{5} \times 4+\frac{1}{5} \times 5+\frac{1}{5} \times 6
|
||||
\\\\
|
||||
&=3.8
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
观察式 1，在定义状态价值函数时，数学期望对于 $G_t$ 没有定义权重或概率，所以每一幕的 $G_t$ 值都是同等价值的，因此，状态价值函数就是很多幕的 $G_t$ 的算术平均值。
|
||||
|
||||
以表 1 中的数据为例：$V(C1)=[3 +(-9)+0+(-13)]/4=-4.75$
|
||||
|
||||
但是，只有 4 幕采样并不能准确计算出真正的期望值，因此，我们需要更多的采样。一般情况下，采样的数量级应该是成千上万的，才会得到一个比较稳定的数学期望值。
|
||||
|
||||
### 4 搭建模型环境
|
||||
|
||||
本部分的代码在 StudentDataModel.py 中。
|
||||
|
||||
#### 定义状态集
|
||||
|
||||
```Python
|
||||
# 状态
|
||||
class States(Enum):
|
||||
Game = 0
|
||||
Class1 = 1
|
||||
Class2 = 2
|
||||
Class3 = 3
|
||||
Pass = 4
|
||||
Rest = 5
|
||||
End = 6
|
||||
```
|
||||
由于本问题中状态比较少，所以可以用枚举方式来定义状态集。
|
||||
|
||||
#### 定义奖励函数
|
||||
|
||||
```Python
|
||||
# 奖励向量
|
||||
# [Game, Class1, Class2, Class3, Pass, Rest, End]
|
||||
Rewards = [-1, -2, -2, -2, 10, 1, 0]
|
||||
```
|
||||
|
||||
#### 定义状态转移矩阵
|
||||
|
||||
```Python
|
||||
# 状态转移概率 from->to
|
||||
P = np.array(
|
||||
[ #Game C1 C2 C3 Pass Rest End
|
||||
[0.9, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.2],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.6, 0.4, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.2, 0.4, 0.4, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
|
||||
]
|
||||
)
|
||||
```
|
||||
这是完全按照图 1 中的标记来定义的，请读者自己核对一下。首先要确定每行的数值的和为 1.0，其次要确定 "from->to" 坐标位置是否正确。万一搞错的话，会给后面写代码时 debug 带来困难。
|
||||
|
||||
在此使用一个简单的向量来定义奖励值，按顺序对应到状态上。
|
||||
|
||||
|
||||
|
||||
#### 定义模型
|
||||
|
||||
按理说有了上面的状态集、奖励、转移矩阵，就可以开始计算 G 值了，但是定义一个统一的模型，会让代码可读性好，出错概率低，并帮助读者加深对概念的理解。
|
||||
|
||||
```Python
|
||||
class DataModel(object):
|
||||
def __init__(self):
|
||||
self.P = P # 状态转移矩阵
|
||||
self.R = Rewards # 奖励
|
||||
self.S = States # 状态集
|
||||
self.num_states = len(self.S) # 状态数量
|
||||
self.end_states = [self.S.End] # 终止状态集
|
||||
|
||||
# 判断给定状态是否为终止状态
|
||||
def is_end(self, s):
|
||||
if (s in self.end_states):
|
||||
return True
|
||||
return False
|
||||
|
||||
# 获得即时奖励，保留此函数可以为将来更复杂的奖励函数做准备
|
||||
def get_reward(self, s):
|
||||
return self.R[s.value]
|
||||
|
||||
# 根据转移概率前进一步，返回（下一个状态、即时奖励、是否为终止）
|
||||
def step(self, curr_s):
|
||||
next_s = np.random.choice(self.S, p=self.P[curr_s.value])
|
||||
return next_s, self.get_reward(next_s), self.is_end(next_s)
|
||||
```
|
||||
|
||||
上面的代码中的注释已经足够丰富了，不再赘述。唯一要提醒的是，我们使用了枚举定义状态，在函数之间传值时都用枚举变量而非具体数值。在函数内部要注意使用 s.value 来做具体索引值。
|
||||
|
||||
当然，如果把奖励函数定义为一个字典，可以直接使用 Reward[State] 的方式来获得当前奖励值，更具可读性。
|
||||
|
||||
### 5 计算状态价值函数
|
||||
|
||||
本部分的代码在 Sampling.py 中。
|
||||
|
||||
OK! 在上一小节，我们的模型已经建立好了，现在可以开始根据式 1 来计算学生学习模型的状态价值函数了。
|
||||
|
||||
#### 算法伪代码
|
||||
|
||||
----
|
||||
输入：起始状态 $S, Episodes, \gamma$
|
||||
多幕 $Episodes$ 循环：
|
||||
$G_{mean} = 0$
|
||||
获得状态 $S$ 的奖励值，看作是 $R_{t+1}$
|
||||
$G \leftarrow R_{t+1} $
|
||||
计数器 $t=1$
|
||||
幕内循环直到终止状态：
|
||||
从 $S$ 根据状态转移概率得到 $S', R'$
|
||||
$G \leftarrow G + \gamma^t R'$
|
||||
$t \leftarrow t+1$
|
||||
$S \leftarrow S'$
|
||||
$G_{mean} \leftarrow G_{mean}+G$
|
||||
返回 $G_{mean} / Episodes$
|
||||
|
||||
---
|
||||
|
||||
#### 算法实现
|
||||
|
||||
```Python
|
||||
def Sampling(dataModel, start_state, episodes, gamma):
|
||||
G_mean = 0 # 定义最终的返回值，G 的平均数
|
||||
# 循环多幕
|
||||
for episode in tqdm.trange(episodes):
|
||||
curr_s = start_state # 把给定的起始状态作为当前状态
|
||||
G = dataModel.get_reward(curr_s) # 由于使用了注重结果奖励方式，所以起始状态也有奖励
|
||||
t = 1 # 折扣因子
|
||||
done = False # 分幕结束标志
|
||||
while (done is False): # 本幕循环
|
||||
next_s, r, done = dataModel.step(curr_s) # 根据当前状态和转移概率获得下一个状态及奖励
|
||||
G += math.pow(gamma, t) * r
|
||||
t += 1
|
||||
curr_s = next_s
|
||||
# end while
|
||||
G_mean += G # 先暂时不计算平均值，而是简单地累加
|
||||
# end for
|
||||
v = G_mean / episodes # 最后再一次性计算平均值，避免增加计算开销
|
||||
return v
|
||||
```
|
||||
|
||||
上述代码可以通过多次循环（由 Episodes）指定，计算指定状态 start_state 的多个回报值 $G$ 的平均值，作为理论上的数学期望值。
|
||||
|
||||
那么 Episodes 的具体数值是多少合适呢？
|
||||
|
||||
#### 多进程并发计算
|
||||
|
||||
在本问题的状态集中一共有 7 个状态。根据上面的算法，首先要指定起始状态，可以遍历状态集中的每个状态作为起始状态。在计算两个状态的状态函数值时互相不干扰，所以，可以考虑使用多进程来并发计算每个指定的起始状态。
|
||||
|
||||
```Python
|
||||
def Sampling_MultiProcess(dataModel, episodes, gamma):
|
||||
pool = mp.Pool(processes=4) # 指定合适的进程数量
|
||||
V = np.zeros((dataModel.num_states))
|
||||
results = []
|
||||
for start_state in dataModel.S: # 遍历状态集中的每个状态作为起始状态
|
||||
results.append(pool.apply_async(Sampling,
|
||||
args=(dataModel, start_state, episodes, gamma,)
|
||||
)
|
||||
)
|
||||
pool.close()
|
||||
pool.join()
|
||||
for s in range(dataModel.num_states):
|
||||
v = results[s].get()
|
||||
V[s] = v
|
||||
|
||||
return V
|
||||
```
|
||||
读者可以根据自己的计算机的 CPU 数量修改 processes=4 的值，但是一定要注意，这个值如果大于你的计算机的 CPU 数量，程序运行速度反而会变慢，因为要在进程间不断切换。
|
||||
|
||||
#### 主过程调用
|
||||
|
||||
```Python
|
||||
if __name__=="__main__":
|
||||
episodes = 10000 # 计算 10000 次的试验的均值作为数学期望值
|
||||
gammas = [0, 0.9, 1] # 指定多个折扣因子做试验
|
||||
dataModel = data.Model()
|
||||
for gamma in gammas:
|
||||
V = Sampling_MultiProcess(dataModel, episodes, gamma)
|
||||
print("gamma =", gamma)
|
||||
for s in dataModel.S:
|
||||
print(str.format("{0}:\t{1}", s.name, V[s.value]))
|
||||
```
|
||||
#### 计算结果
|
||||
```
|
||||
gamma = 0
|
||||
Game: -1.0
|
||||
Class1: -2.0
|
||||
Class2: -2.0
|
||||
Class3: -2.0
|
||||
Pass: 10.0
|
||||
Rest: 1.0
|
||||
End: 0.0
|
||||
---------------
|
||||
gamma = 0.9
|
||||
Game: -7.619960436758965
|
||||
Class1: -5.012940670268943
|
||||
Class2: 0.9443617132189657
|
||||
Class3: 4.020495688121447
|
||||
Pass: 10.0
|
||||
Rest: 1.908911196549832
|
||||
End: 0.0
|
||||
---------------
|
||||
gamma = 1
|
||||
Game: -22.1306
|
||||
Class1: -12.638
|
||||
Class2: 1.3951
|
||||
Class3: 4.4477
|
||||
Pass: 10.0
|
||||
Rest: 0.8086
|
||||
End: 0.0
|
||||
```
|
||||
|
||||
数据解读：
|
||||
|
||||
- $\gamma=0$ 时
|
||||
|
||||
价值函数值就等于当前状态的奖励值，这和价值函数以及回报值的定义相符。
|
||||
|
||||
- $\gamma=1$ 时
|
||||
|
||||
距离终点越远的状态（比如 Game 和 Class1），其价值函数越小，甚至都不在一个数量级上了。
|
||||
|
||||
- $\gamma=0.9$ 时
|
||||
|
||||
由于折扣的存在，使得各个状态的价值函数值之间有些许的平衡。
|
||||
|
||||
### 6 如何确定采样的次数
|
||||
|
||||
具体地说就是如何确定算法中的幕数 Episodes 的数值。
|
||||
|
||||
根据问题的复杂程度不同，幕数必然会不同。但到目前为止，没有人从理论层面研究过这个问题，所以有一些偏实践的方法，供大家参考。
|
||||
|
||||
#### 试探
|
||||
|
||||
先用比较小的数值，比如 100，去做几次尝试，如果发现几次尝试的结果之间有很大的方差，就增加到 1000 再试试 ...... 以此类推，也许到 10 万时才能相对稳定。
|
||||
|
||||
方差公式为：
|
||||
|
||||
$$
|
||||
\sigma^2=\frac{1}{n}\sum_{i=1}^n (V_i-\mu)^2 \tag{2}
|
||||
$$
|
||||
|
||||
比如，我们只关注 Rest 状态的价值函数值，运行 3 次得到 $V_1,V_2,V_3$ 的值，而 $\mu=\frac{1}{3}(V_1+V_2+V_3)$ 是均值：
|
||||
|
||||
1. 当 episodes=100 时，运行 3 次，结果分别是 [1.2, 2.8, 2.0]，相差非常大，$\sigma^2=0.427$；
|
||||
2. 设置 episodes=1000，运行 3 次，结果分别是 [0.5, 1.2, 0.9]，方差减小了，$\sigma^2=0.082$，但还不够好；
|
||||
3. 设置 episodes=10000，运行 3 次，结果分别是 [0.86, 0.75, 0.82]，$\sigma^2=0.002$，如果这个方差到了你的心理预期，就可以结束了，否则可以再增加 episodes 的次数。
|
||||
|
||||
#### 比较
|
||||
|
||||
先看一个增量计算平均值的公式：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
V_{n+1} &= \frac{1}{n+1} \sum_{i=1}^{n+1} G_i
|
||||
\\\\
|
||||
&= \frac{1}{n+1}(G_{n+1}+nV_n)= \frac{1}{n+1}(G_{n+1}+n V_n+ V_n-V_n)
|
||||
\\\\
|
||||
&=V_n + \frac{1}{n+1}(G_{n+1}-V_n)
|
||||
\end{aligned}
|
||||
\tag{3}
|
||||
$$
|
||||
|
||||
式 3 表达的意思是，n+1 幕时 $G_{n+1}$ 的期望值 $V_{n+1}$，等于 $n$ 幕时的 $G_n$ 期望值 $V_{n}$，再加上 n+1 幕时的 $G_{n+1}$ 与 $V_n$ 的差值除以 (n+1)。
|
||||
|
||||
把式 2 变形得到：
|
||||
|
||||
$$
|
||||
V_{n+1} - V_{n} = \frac{1}{n+1}(G_{n+1}-V_n)
|
||||
\tag{4}
|
||||
$$
|
||||
|
||||
|
||||
在第 n+1 幕时，先计算出式 4 的等号后面的部分，检查这个值是否足够小（比如小于 1e-2），就可以认为已经收敛了。
|
||||
|
||||
当然，这里的步长可以不是 1 幕，而是 10 幕或者 100 幕。
|
||||
|
||||
### 问题与讨论
|
||||
|
||||
1. 为什么在不同的折扣因子情况下，Pass 状态的值永远是 10.0，而 End 状态的值永远是 0.0 ？
|
||||
2. 请使用试探法来找到比较理想的分幕次数。
|
||||
3. 请使用比较法来找到比较理想的分幕次数。
|
||||
4. 假设 $Episodes=10000,\gamma=0.9$，多次计算状态值，哪个状态的方差最大？为什么？
|
||||
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 41 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 26 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 52 KiB |
@@ -0,0 +1,79 @@
|
||||
|
||||
from asyncio import create_subprocess_shell
|
||||
import numpy as np
|
||||
from enum import Enum
|
||||
|
||||
# 状态
|
||||
class States(Enum):
|
||||
Start = 0 # 出发
|
||||
Normal = 1 # 正常行驶
|
||||
Pedestrians = 2 # 礼让行人
|
||||
DownSpeed = 3 # 闹市减速
|
||||
ExceedSpeed = 4 # 超速行驶
|
||||
RedLight = 5 # 路口闯灯
|
||||
LowSpeed = 6 # 小区减速
|
||||
MobilePhone = 7 # 拨打手机
|
||||
Crash = 8 # 发生事故
|
||||
Goal = 9 # 安全抵达
|
||||
End = 10 # 结束
|
||||
|
||||
# 状态转移概率
|
||||
P = np.array(
|
||||
[
|
||||
[0.0, 0.9, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.2, 0.1, 0.1, 0.1, 0.3, 0.1, 0.0, 0.1, 0.0],
|
||||
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.7, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.2, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.5, 0.0, 0.0],
|
||||
[0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
|
||||
[0.0, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
|
||||
]
|
||||
)
|
||||
|
||||
# 奖励向量
|
||||
# |出发|正常行驶|礼让行人|闹市减速|超速行驶|路口闯灯|小区减速|拨打手机|发生事故|安全抵达|结束|
|
||||
R = [0, 0, +1, +1, -3, -6, +1, -3, -1, +5, 0]
|
||||
|
||||
class DataModel(object):
|
||||
def __init__(self):
|
||||
self.P = P # 状态转移矩阵
|
||||
self.R = R # 奖励
|
||||
self.S = States # 状态集
|
||||
self.num_states = len(self.S) # 状态数量
|
||||
self.end_states = [self.S.End] # 终止状态集
|
||||
self.V_ground_truth = Matrix(self, 1)
|
||||
|
||||
# 判断给定状态是否为终止状态
|
||||
def is_end(self, s):
|
||||
if (s in self.end_states):
|
||||
return True
|
||||
return False
|
||||
|
||||
# 获得即时奖励，保留此函数可以为将来更复杂的奖励函数做准备
|
||||
def get_reward(self, s):
|
||||
return self.R[s.value]
|
||||
|
||||
# 根据转移概率前进一步，返回（下一个状态、即时奖励、是否为终止）
|
||||
def step(self, curr_s):
|
||||
next_s = np.random.choice(self.S, p=self.P[curr_s.value])
|
||||
return next_s, self.get_reward(next_s), self.is_end(next_s)
|
||||
|
||||
def Matrix(dataModel, gamma):
|
||||
num_state = dataModel.P.shape[0]
|
||||
I = np.eye(dataModel.num_states) * (1+1e-7)
|
||||
#I = np.eye(dataModel.num_states)
|
||||
tmp1 = I - gamma * dataModel.P
|
||||
tmp2 = np.linalg.inv(tmp1)
|
||||
vs = np.dot(tmp2, dataModel.R)
|
||||
return vs
|
||||
|
||||
|
||||
|
||||
if __name__=="__main__":
|
||||
dataModel = DataModel()
|
||||
v = Matrix(dataModel, 1.0)
|
||||
print(np.around(v,2))
|
||||
+18
-10
@@ -3,7 +3,7 @@ import tqdm
|
||||
import multiprocessing as mp
|
||||
import math
|
||||
import numpy as np
|
||||
import StudentDataModel as data
|
||||
import DriveDataModel as data
|
||||
import time
|
||||
|
||||
def Sampling_MultiProcess(dataModel, episodes, gamma):
|
||||
@@ -23,36 +23,44 @@ def Sampling_MultiProcess(dataModel, episodes, gamma):
|
||||
|
||||
return V
|
||||
|
||||
|
||||
# 多次采样获得回报 G 的数学期望，即状态价值函数 V
|
||||
def Sampling(dataModel, start_state, episodes, gamma):
|
||||
G_mean = 0 # 定义最终的返回值，G 的平均数
|
||||
G_sum = 0 # 定义最终的返回值，G 的平均数
|
||||
# 循环多幕
|
||||
for episode in tqdm.trange(episodes):
|
||||
# 由于使用了注重结果奖励方式，所以起始状态也有奖励，做为 G 的初始值
|
||||
G = dataModel.get_reward(start_state)
|
||||
curr_s = start_state # 把给定的起始状态作为当前状态
|
||||
G = dataModel.get_reward(curr_s) # 由于使用了注重结果奖励方式，所以起始状态也有奖励
|
||||
t = 1 # 折扣因子
|
||||
t = 1 # 折扣因子
|
||||
done = False # 分幕结束标志
|
||||
while (done is False): # 本幕循环
|
||||
next_s, r, done = dataModel.step(curr_s) # 根据当前状态和转移概率获得下一个状态及奖励
|
||||
# 根据当前状态和转移概率获得:下一个状态,奖励,是否到达终止状态
|
||||
next_s, r, done = dataModel.step(curr_s)
|
||||
G += math.pow(gamma, t) * r
|
||||
t += 1
|
||||
curr_s = next_s
|
||||
# end while
|
||||
G_mean += G # 先暂时不计算平均值，而是简单地累加
|
||||
G_sum += G # 先暂时不计算平均值，而是简单地累加
|
||||
# end for
|
||||
v = G_mean / episodes # 最后再一次性计算平均值，避免增加计算开销
|
||||
return v
|
||||
V = G_sum / episodes # 最后再一次性计算平均值，避免增加计算开销
|
||||
return V
|
||||
|
||||
def RMSE(a,b):
|
||||
err = np.sqrt(np.sum(np.square(a - b))/a.shape[0])
|
||||
return err
|
||||
|
||||
if __name__=="__main__":
|
||||
start = time.time()
|
||||
episodes = 10000 # 计算 10000 次的试验的均值作为数学期望值
|
||||
gammas = [0, 0.9, 1] # 指定多个折扣因子做试验
|
||||
Vs = []
|
||||
dataModel = data.DataModel()
|
||||
for gamma in gammas:
|
||||
V = Sampling_MultiProcess(dataModel, episodes, gamma)
|
||||
Vs.append(V)
|
||||
print("gamma =", gamma)
|
||||
for s in dataModel.S:
|
||||
print(str.format("{0}:\t{1}", s.name, V[s.value]))
|
||||
end = time.time()
|
||||
print(end-start)
|
||||
#print(end-start)
|
||||
print("RMSE = ", RMSE(Vs[2], dataModel.V_ground_truth))
|
||||
@@ -0,0 +1,86 @@
|
||||
|
||||
import tqdm
|
||||
import multiprocessing as mp
|
||||
import math
|
||||
import numpy as np
|
||||
import DriveDataModel as data
|
||||
import time
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
def Sampling_MultiProcess(dataModel, episodes, gamma, checkpoint):
|
||||
pool = mp.Pool(processes=4) # 指定合适的进程数量
|
||||
V = dict()
|
||||
results = []
|
||||
for start_state in dataModel.S: # 遍历状态集中的每个状态作为起始状态
|
||||
results.append(pool.apply_async(Sampling_Checkpoint,
|
||||
args=(dataModel, start_state, episodes, gamma, checkpoint,)
|
||||
)
|
||||
)
|
||||
pool.close()
|
||||
pool.join()
|
||||
for s in range(dataModel.num_states):
|
||||
v = results[s].get()
|
||||
V[s] = v
|
||||
|
||||
return V
|
||||
|
||||
# 多次采样获得回报 G 的数学期望，即状态价值函数 V
|
||||
def Sampling_Checkpoint(dataModel, start_state, episodes, gamma, checkpoint):
|
||||
V = []
|
||||
G_sum = 0 # 定义最终的返回值，G 的平均数
|
||||
# 循环多幕
|
||||
for episode in tqdm.trange(episodes):
|
||||
# 由于使用了注重结果奖励方式，所以起始状态也有奖励，做为 G 的初始值
|
||||
G = dataModel.get_reward(start_state)
|
||||
curr_s = start_state # 把给定的起始状态作为当前状态
|
||||
t = 1 # 折扣因子
|
||||
done = False # 分幕结束标志
|
||||
while (done is False): # 本幕循环
|
||||
# 根据当前状态和转移概率获得:下一个状态,奖励,是否到达终止状态
|
||||
next_s, r, done = dataModel.step(curr_s)
|
||||
G += math.pow(gamma, t) * r
|
||||
t += 1
|
||||
curr_s = next_s
|
||||
# end while
|
||||
G_sum += G # 先暂时不计算平均值，而是简单地累加
|
||||
if (episode+1)%checkpoint == 0:
|
||||
V.append(G_sum / (episode+1))
|
||||
# end for
|
||||
V.append(G_sum / episodes)
|
||||
return V
|
||||
|
||||
|
||||
def RMSE(a,b):
|
||||
err = np.sqrt(np.sum(np.square(a - b))/a.shape[0])
|
||||
return err
|
||||
|
||||
def test_once():
|
||||
episodes = 10000 # 计算 10000 次的试验的均值作为数学期望值
|
||||
gamma = 1
|
||||
checkpoint = 100
|
||||
dataModel = data.DataModel()
|
||||
V = Sampling_MultiProcess(dataModel, episodes, gamma, checkpoint)
|
||||
num_checkpoint = len(V[0])
|
||||
array = np.zeros((num_checkpoint, dataModel.num_states))
|
||||
for s_value in range(dataModel.num_states):
|
||||
array[:,s_value] = V[s_value] # V is dictionary
|
||||
|
||||
errors = []
|
||||
for i in range(num_checkpoint):
|
||||
err = RMSE(array[i], dataModel.V_ground_truth)
|
||||
errors.append(err)
|
||||
|
||||
return errors
|
||||
|
||||
if __name__=="__main__":
|
||||
ERRORS = []
|
||||
for i in range(10):
|
||||
errors = test_once()
|
||||
ERRORS.append(errors)
|
||||
|
||||
avg_E = np.mean(ERRORS, axis=0)
|
||||
plt.plot(avg_E)
|
||||
plt.grid()
|
||||
plt.title(str.format("min RMSE={0}", np.min(avg_E)))
|
||||
plt.show()
|
||||
|
||||
@@ -0,0 +1,336 @@
|
||||
## 状态价值函数
|
||||
|
||||
### 1 提出问题
|
||||
|
||||
在安全驾驶问题中，我们学习了马尔可夫奖励过程，给其中的各个状态以奖励值，或正或负。比如：
|
||||
|
||||
- 礼让行人、闹市减速、小区减速等状态都可以得到 +1 的奖励。
|
||||
- 开车打电话、闯红灯、超速行驶等状态可以得到 -3 或 -6 的扣分。
|
||||
|
||||
读者不禁会产生几个问题：
|
||||
|
||||
1. 这个奖励值可以真正表示这个状态的好坏吗？
|
||||
2. 具有相同奖励值的状态，哪个更好？
|
||||
3. 虽然有回报值 $G$ 的定义，但是由于分幕数据序列的不同，对某个状态来说，每次采样得到的 $G$ 值都不一样，以哪一个为准呢？
|
||||
|
||||
### 2 建立模型
|
||||
|
||||
在学习马尔可夫奖励过程时，已经建立好了这个问题的状态转移模型和奖励模型，可以直接拿过来用。
|
||||
|
||||
图 1 是该问题的马尔可夫奖励模型，每个状态都有一个即时奖励值。
|
||||
|
||||
<center>
|
||||
<img src="./img/Drive-2.png" width="500">
|
||||
|
||||
图 1 学习问题的状态转移概率图
|
||||
</center>
|
||||
|
||||
需要再次说明的是，在图 1 中，我们使用了**注重结果**的奖励定义方式，直接给每个状态赋值一个奖励，意味只要达到这个状态，就可以立刻获得标注出的奖励值，而不管是从哪条路径达到的。
|
||||
|
||||
表 1 中列出了状态转移矩阵，与租车问题中的矩阵形式相同。
|
||||
|
||||
表 1 状态转移矩阵
|
||||
|
||||
|P: 从$\rightarrow$到|出发|正常<br>行驶|礼让<br>行人|闹市<br>减速|超速<br>行驶|路口<br>闯灯|小区<br>减速|拨打<br>手机|发生<br>事故|安全<br>抵达|结束|
|
||||
|-|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|
||||
|S：出发|-|0.9|-|-|0.1|-|-|-|-|-|-|
|
||||
|N：正常行驶|-|-|0.2|0.1|0.1|0.1|0.3|0.1|-|0.1|-|
|
||||
|P：礼让行人|-|1.0|-|-|-|-|-|-|-|-|-|
|
||||
|D：闹市减速|-|0.7|0.3|-|-|-|-|-|-|-|-|
|
||||
|X：超速行驶|-|0.2|-|-|-|0.3|-|-|0.5|-|-|
|
||||
|R：路口闯灯|-|0.1|-|-|-|-|-|-|0.9|-|-|
|
||||
|L：小区减速|-|-|-|-|-|-|-|-|-|1.0|-|
|
||||
|M：拨打手机|-|0.6|-|-|-|-|-|-|0.4|-|-|
|
||||
|C：发生事故|-|-|-|-|-|-|-|-|-|-|1.0|
|
||||
|G：安全抵达|-|-|-|-|-|-|-|-|-|-|1.0|
|
||||
|E：结束|-|-|-|-|-|-|-|-|-|-|1.0|
|
||||
|
||||
|
||||
有的读者可能会较真儿：进入小区减速状态后，就一定可以安全抵达吗？当然，在小区里可能会有很多突发情况，比如老人儿童突然横穿道路等等。在此我们就不再细化这个模型了。
|
||||
|
||||
在上一节中，已经有了分幕、奖励、回报的概念，这一节中，将会利用这些基础概念来定义每个**状态价值函数**，从而可以比较状态之间的好坏。
|
||||
|
||||
### 3 状态价值函数（State Value Function）
|
||||
|
||||
在上一节通过对 $G$ 的计算，以及对状态图的分析理解，我们似乎已经得到了一些启示：距离终点越近的状态，越接近于成功，它的状态价值就越高，似乎用 $G$ 值就可以表示该时刻的状态价值。
|
||||
|
||||
但是，会有两个麻烦出现：
|
||||
|
||||
1. 每个司机经过的路径都不完全一样，状态图虽然是有向的，但是由于环状转移的原因，一个状态在一个序列中可能会出现多次。比如 *正常行驶* 状态，有 6 种情况可以到达；
|
||||
2. 从 *出发* 开始的路径又可以有很多种，如表 2 所示。
|
||||
|
||||
表 2 以 *出发* 开始的 $G$ 值的计算
|
||||
|
||||
||分幕采样序列|回报值计算（$\gamma=1$）|
|
||||
|-|-|-|
|
||||
|1|Start - N - L - G - End|$G_{S}=0+0+1+5+0=6$|
|
||||
|2|Start - N - P - N - L - G - End|$G_{S}=0+0+1+0+1+5+0=7$|
|
||||
|3|Start - N - R - C - End|$G_{S}=0+0-6-1+0=-7$|
|
||||
|4|Start - X - R - N - R - C - End|$G_{S}=0-3-6+0-6-1+0=-16$|
|
||||
|
||||
|
||||
这样的话，一个状态就可能有 4 个不同的 $G$ 值，我们用哪个当作其价值函数呢？另外，仔细观察 $G$ 的表达式，它只与时刻及奖励有关，没有体现出状态来。
|
||||
|
||||
考虑到以上两点，定义状态价值函数如下：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
V_t(s) &= \mathbb E [G_t | S_t = s]
|
||||
\\\\
|
||||
&=\mathbb E [ R_{t+1}+\gamma R_{t+2}+\gamma^2 R_{t+3}+ \gamma^3 R_{t+4}+ \cdots]
|
||||
\end{aligned}
|
||||
\tag{1}
|
||||
$$
|
||||
|
||||
式 1 的含义是，定义状态价值函数 $V$ 是回报 $G$ 的**数学期望**。时刻 $t$ 在这里只起到一个按顺序串联状态 $S$，从而得到奖励 $R$ 的作用。
|
||||
|
||||
### 4 数学期望
|
||||
|
||||
简单地回忆一下**数学期望**的概念。
|
||||
|
||||
一个正常的六面的骰子，投出去后可以得到 [1,2,3,4,5,6] 六种结果，而且概率相等，那么这个骰子的期望值是 $(1+2+3+4+5+6)/6=3.5$。哈哈，读者可能会发现 3.5 这个数子，骰子无法投出来，所以它只是一种定义。
|
||||
|
||||
但是，一个不正常的骰子，比如 [4,5,6] 出现的概率 $p$ 都是 $\frac{1}{5}$，而 [1,2,3] 出现的概率 $p$ 都是 $\frac{2}{15}$，那么它的数学期望是：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
\mathbb E[骰子]&=\sum_{i=1}^6 p_i V_i
|
||||
\\\\
|
||||
&= \frac{2}{15} \times 1+\frac{2}{15} \times 2+\frac{2}{15} \times 3+\frac{1}{5} \times 4+\frac{1}{5} \times 5+\frac{1}{5} \times 6
|
||||
\\\\
|
||||
&=3.8
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
观察式 2，在定义状态价值函数时，数学期望对于 $G_t$ 没有定义权重或概率，所以每一幕的 $G_t$ 值都是同等价值的，因此，状态价值函数就是很多幕的 $G_t$ 的算术平均值。
|
||||
|
||||
以表 1 中的数据为例：$V(Start)=[6+7+(-7)+(-16)]/4=-2.5$
|
||||
|
||||
但是，只有 4 幕采样并不能准确计算出真正的期望值，因此，我们需要更多的采样。一般情况下，采样的数量级应该是成千上万的，才会得到一个比较稳定的数学期望值。
|
||||
|
||||
### 5 搭建模型环境
|
||||
|
||||
【代码位置：DriveDataModel.py】
|
||||
|
||||
#### 定义状态集
|
||||
|
||||
```Python
|
||||
# 状态
|
||||
class States(Enum):
|
||||
Start = 0 # 出发
|
||||
Normal = 1 # 正常行驶
|
||||
Pedestrians = 2 # 礼让行人
|
||||
DownSpeed = 3 # 闹市减速
|
||||
ExceedSpeed = 4 # 超速行驶
|
||||
RedLight = 5 # 路口闯灯
|
||||
LowSpeed = 6 # 小区减速
|
||||
MobilePhone = 7 # 拨打手机
|
||||
Crash = 8 # 发生事故
|
||||
Goal = 9 # 安全抵达
|
||||
End = 10 # 结束
|
||||
```
|
||||
由于本问题中状态比较少，所以可以用枚举方式来定义状态集。
|
||||
|
||||
#### 定义奖励函数
|
||||
|
||||
```Python
|
||||
# 奖励向量
|
||||
# |出发|正常行驶|礼让行人|闹市减速|超速行驶|路口闯灯|小区减速|拨打手机|发生事故|安全抵达|结束|
|
||||
R = [0, 0, +1, +1, -3, -6, +1, -3, -1, +5, 0]
|
||||
```
|
||||
|
||||
#### 定义状态转移矩阵
|
||||
|
||||
```Python
|
||||
# 状态转移概率 from->to
|
||||
P = np.array(
|
||||
[
|
||||
[0.0, 0.9, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.2, 0.1, 0.1, 0.1, 0.3, 0.1, 0.0, 0.1, 0.0],
|
||||
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.7, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.2, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.5, 0.0, 0.0],
|
||||
[0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
|
||||
[0.0, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
|
||||
]
|
||||
)
|
||||
```
|
||||
这是完全按照图 1 中的标记来定义的，请读者自己核对一下。首先要确定每行的数值的和为 1.0，其次要确定 "from->to" 坐标位置是否正确。万一搞错的话，会给后面写代码时 debug 带来困难。
|
||||
|
||||
在此使用一个简单的向量来定义奖励值，按顺序对应到状态上。
|
||||
|
||||
|
||||
|
||||
#### 定义模型
|
||||
|
||||
按理说有了上面的状态集、奖励、转移矩阵，就可以开始计算 G 值了。但是定义一个统一的模型，把细节封装成一些标准的交互函数，会让代码可读性好，出错概率低，并帮助读者加深对概念的理解。
|
||||
|
||||
```Python
|
||||
class DataModel(object):
|
||||
def __init__(self):
|
||||
self.P = P # 状态转移矩阵
|
||||
self.R = R # 奖励
|
||||
self.S = States # 状态集
|
||||
self.num_states = len(self.S) # 状态数量
|
||||
self.end_states = [self.S.End] # 终止状态集
|
||||
|
||||
# 判断给定状态是否为终止状态
|
||||
def is_end(self, s):
|
||||
if (s in self.end_states):
|
||||
return True
|
||||
return False
|
||||
|
||||
# 获得即时奖励，保留此函数可以为将来更复杂的奖励函数做准备
|
||||
def get_reward(self, s):
|
||||
return self.R[s.value]
|
||||
|
||||
# 根据转移概率前进一步，返回（下一个状态、即时奖励、是否为终止）
|
||||
def step(self, curr_s):
|
||||
next_s = np.random.choice(self.S, p=self.P[curr_s.value])
|
||||
return next_s, self.get_reward(next_s), self.is_end(next_s)
|
||||
```
|
||||
|
||||
上面的代码中的注释已经足够丰富了，不再赘述。唯一要提醒的是，我们使用了枚举定义状态，在函数之间传值时都用枚举变量而非具体数值。在函数内部要注意使用 s.value 来做具体索引值。
|
||||
|
||||
当然，如果把奖励函数定义为一个字典，可以直接使用 Reward[State] 的方式来获得当前奖励值，更具可读性。
|
||||
|
||||
### 6 计算状态价值函数
|
||||
|
||||
【代码位置：DriveSampling.py】
|
||||
|
||||
OK! 模型已经建立好了，现在可以开始根据式 1 来计算学生学习模型的状态价值函数了。
|
||||
|
||||
#### 算法描述
|
||||
|
||||
----
|
||||
输入：起始状态 $S, Episodes, \gamma$
|
||||
初始化：$G_{sum} \leftarrow 0$ # 累计多幕的G值以便求平均
|
||||
多幕 $Episodes$ 循环：
|
||||
获得状态 $S$ 的奖励值，看作是 $R_{t+1}$
|
||||
$G \leftarrow R_{t+1} $
|
||||
计数器 $t=1$
|
||||
幕内循环直到终止状态：
|
||||
从 $S$ 根据状态转移概率得到 $S',R'$ 以及终止标志
|
||||
$G \leftarrow G + \gamma^t R'$
|
||||
$t \leftarrow t+1$
|
||||
$S \leftarrow S'$
|
||||
$G_{sum} \leftarrow G_{sum}+G$
|
||||
$V_S \leftarrow G_{sum} / Episodes$
|
||||
输出：$V_S$
|
||||
|
||||
---
|
||||
|
||||
#### 算法说明
|
||||
|
||||
以表 2 中的第 1 个采样序列为例，说明算法的执行过程。
|
||||
|
||||
<center>
|
||||
<img src="./img/Sampling.png" width="80%">
|
||||
|
||||
图 2 采样算法说明
|
||||
</center>
|
||||
|
||||
图 2 中展示了以5个状态组成的序列为例的 $G$ 的计算过程：
|
||||
|
||||
1. 起始状态为 S，得到奖励 $R_1$，保存到 $G$ 中，其中的下标 $G_{[1]}$ 表示第一步；
|
||||
2. 转移到状态 N，得到奖励 $R_2$，乘以 $\gamma$ 后与第一步的 $G$ 相加，仍然保存到 $G$ 中，原来的 $G$ 值就被替换掉了；
|
||||
3. 以此类推，一直到最后的 E 状态，得到 $R_T$，与第 $[4]$ 步的 $G$ 相加，终止幕内循环，得到状态 S 的一个采样序列的回报值 $G_{S}$；
|
||||
4. 多次重复上述过程，得到不同的采样序列的 $G$ 值，累计；
|
||||
5. 最后的累计值除以幕数，就可以认为是 $G$ 的数学期望，因而得到状态值 $V_S$。
|
||||
|
||||
这个算法的特点是使用了最少的变量，在算法过程中，一共只用了 $G_{sum},G,R,t,S,S'$ 等几个变量，没有使用任何列表或数组，而且完全是按照回报的定义以及状态价值函数的定义来实现的，便于读者理解。
|
||||
|
||||
#### 算法实现
|
||||
|
||||
```Python
|
||||
# 多次采样获得回报 G 的数学期望，即状态价值函数 V
|
||||
def Sampling(dataModel, start_state, episodes, gamma):
|
||||
G_sum = 0 # 定义最终的返回值，G 的平均数
|
||||
# 循环多幕
|
||||
for episode in tqdm.trange(episodes):
|
||||
# 由于使用了注重结果奖励方式，所以起始状态也有奖励，做为 G 的初始值
|
||||
G = dataModel.get_reward(start_state)
|
||||
curr_s = start_state # 把给定的起始状态作为当前状态
|
||||
t = 1 # 折扣因子
|
||||
done = False # 分幕结束标志
|
||||
while (done is False): # 本幕循环
|
||||
# 根据当前状态和转移概率获得:下一个状态,奖励,是否到达终止状态
|
||||
next_s, r, done = dataModel.step(curr_s)
|
||||
G += math.pow(gamma, t) * r
|
||||
t += 1
|
||||
curr_s = next_s
|
||||
# end while
|
||||
G_sum += G # 先暂时不计算平均值，而是简单地累加
|
||||
# end for
|
||||
Vs = G_sum / episodes # 最后再一次性计算平均值，避免增加计算开销
|
||||
return Vs
|
||||
```
|
||||
|
||||
上述代码可以通过多次循环（由 Episodes）指定，计算指定状态 start_state 的多个回报值 $G$ 的平均值，作为理论上的数学期望值。
|
||||
|
||||
那么 Episodes 的具体数值是多少合适呢？后面再解释。
|
||||
|
||||
#### 多进程并发计算
|
||||
|
||||
在本问题的状态集中一共有 11 个状态。根据上面的算法，首先要指定起始状态，可以遍历状态集中的每个状态作为起始状态。在计算两个状态的状态函数值时互相不干扰，所以，可以考虑使用多进程来并发计算每个指定的起始状态。
|
||||
|
||||
读者可以根据自己的计算机的 CPU 数量修改 processes=4 的值，但是一定要注意，这个值如果大于你的计算机的 CPU 数量，程序运行速度反而会变慢，因为要在进程间不断切换。
|
||||
|
||||
#### 主过程调用
|
||||
|
||||
```Python
|
||||
if __name__=="__main__":
|
||||
episodes = 10000 # 计算 10000 次的试验的均值作为数学期望值
|
||||
gammas = [0, 0.9, 1] # 指定多个折扣因子做试验
|
||||
dataModel = data.Model()
|
||||
for gamma in gammas:
|
||||
V = Sampling_MultiProcess(dataModel, episodes, gamma) # 多进程调用
|
||||
print("gamma =", gamma)
|
||||
for s in dataModel.S:
|
||||
print(str.format("{0}:\t{1}", s.name, V[s.value]))
|
||||
```
|
||||
#### 计算结果
|
||||
|
||||
表 各个状态的奖励值和价值函数
|
||||
|
||||
|状态|R|$\gamma = 0$|$\gamma = 0.9$|$\gamma = 1$|
|
||||
|-|-:|:-:|:-:|:-:|
|
||||
|出发 Start| 0| 0.0 |0.47|1.09|
|
||||
|正常行驶 Normal| 0| 0.0 |1.25|1.64|
|
||||
|礼让行人 Pedestrians| +1| 1.0 |2.11|2.74|
|
||||
|闹市减速 DownSpeed| +1| 1.0 |2.39|2.99|
|
||||
|超速行驶 ExceedSpeed| -3| -3.0 |-5.03|-5.18|
|
||||
|路口闯灯 RedLight| -6| -6.0|-6.73|-6.72|
|
||||
|小区减速 LowSpeed| +1| 1.0 |5.5|6.0|
|
||||
|拨打电话 MobilePhone| -3| -3.0|-2.69|-2.32|
|
||||
|发生事故 Crash| -1| -1.0|-1.0|-1.0|
|
||||
|安全抵达 Goal| +5| 5.0 |5.0|5.0|
|
||||
|终止 End| 0| 0.0 |0.0|0.0|
|
||||
|
||||
数据解读：
|
||||
|
||||
- $\gamma=0$ 时，价值函数值就等于当前状态的奖励值，这和价值函数以及回报值的定义相符。
|
||||
- $\gamma \ne 0$ 时，
|
||||
- 正常行驶状态的奖励值虽然为 0，但是其状态价值是大于 0 的。
|
||||
- 有正的奖励的状态，价值函数值都会比奖励值大，除了 *安全抵达* 之外。
|
||||
- 有负的奖励的状态，价值函数值都会比奖励值小，除了 *发生事故* 之外。
|
||||
- $\gamma = 1$ 的值要比 $\gamma = 0.9$ 的值大一些，因为折扣系数的原因。
|
||||
|
||||
一个有趣的问题，可以帮助读者更好地理解状态价值函数：比较 *闹市减速* 和 *小区减速* 两个状态，同样是 *减速* ，为什么后者的状态价值函数值比前者高呢？
|
||||
|
||||
可以从驾驶的过程中这样理解：
|
||||
|
||||
- 在 *闹市减速* 状态时，虽然司机此时做了正确的选择，使车辆处于合理的行驶状态，但并不一定保证在后续的驾驶过程中，司机同样会做出正确的选择，仍然有可能超速、闯灯等等，以至于发生交通事故。
|
||||
- 在 *小区减速* 状态时，从状态图上看，因为后续会以 100% 的概率进入 *安全抵达* 状态，没有变数，可以得到很高的奖励。所以状态价值函数值会比 *闹市减速* 状态要高。
|
||||
|
||||
|
||||
### 问题与讨论
|
||||
|
||||
1. 为什么在不同的折扣因子情况下，*安全抵达* 的状态的值永远是 10.0，而 *终止* 状态的值永远是 0.0 ？
|
||||
2. 假设 $Episodes=10000,\gamma=0.9$，多次计算状态值，哪个状态的方差最大？为什么？
|
||||
3. 为什么 *安全抵达* 的状态值和奖励值相等，而其它的具有正奖励值的状态值都比奖励值大？
|
||||
|
||||
@@ -0,0 +1,226 @@
|
||||
|
||||
## 蒙特卡洛方法
|
||||
|
||||
本节中将要学习的蒙特卡洛方法，是强化学习的一种重要手段，在后面的学习中，要更充分地讨论蒙特卡洛预测和蒙特卡洛控制，都是以本节的内容为理论基础的。
|
||||
|
||||
其实在“三门问题”和“学生学习问题”中，已经使用了蒙特卡洛的思想来解决问题了，接下来我们系统地了解一下这种方法。
|
||||
|
||||
### 1 提出问题一
|
||||
|
||||
#### 估算圆周率
|
||||
|
||||
假设我们已经知道了圆面积的方程为 $S = \pi r^2$，但是不知道 $\pi$ 的具体数值。能否通过一些简单的方法得到 $\pi$ 值呢？
|
||||
|
||||
请读者先放松一下，看看历史上关于圆周率计算的故事。
|
||||
|
||||
联合国教科文组织在2019年11月26日第四十届大会批准宣布，3月14日为“国际数学日（International Day of Mathematics,简称IDM）”。因为“3.14”是圆周率数值最接近的数字，所以这一天也叫圆周率日（$\pi$ day）。
|
||||
|
||||
- 公元前250年，希腊数学家阿基米德通过割圆术计算圆周率，阿基米德进行了96边形的割圆之后，将圆周率推到了小数点后两位3.14。
|
||||
- 直到公元265年，中国的数学家刘徽用割圆术的方法，通过正3072边形计算出π的数值为3.1416，艰难地把圆周率推到了小数点后四位。
|
||||
- 200年后，祖冲之继续使用割圆术计算12,288形的边长，将圆周率推到了小数点后六位，可惜的是，由于文献的失传，祖冲之的计算方法我们现在已经不得而知了。
|
||||
- 800年后，随着近代数学的发展，数学家韦达、罗门、科伊伦、司乃耳、格林伯格通过割圆术陆续将圆周率推到了小数点后39位，这个精度是什么概念呢，如果我们通过小数点后39位的圆周率计算一个一个可观察宇宙大小的圆，计算的误差仅仅只有一个氢原子大小。
|
||||
- 十六世纪到十七世纪，人们发现了一种新的圆周率计算方法——无穷级数法，让计算圆周率的工作变得更加快速。无穷级数是一组无穷数列的和，数学家梅钦通过无穷级数将圆周率推算到小数点后100位，在很短的时间里，人们通过梅钦类公式反复打破了新的圆周率记录。
|
||||
- 18世纪，法国数学家布丰提出了随机投针法，即利用概率统计的方法来计算圆周率的值，也就是著名的投针实验。布丰在地板上画出若干平行的直线，再将一根根短于平行直线距离的针撒到地板上，通过统计针的总数和与直线相交的针的个数，从而计算圆周率。
|
||||
|
||||
这种算法虽然虽然没有打破圆周率的记录，但这种将几何与概率结合起来的思想催生了蒙特卡洛算法，也让人工智能成为了可能。
|
||||
|
||||
#### 用蒙特卡洛方法估算圆周率
|
||||
|
||||
先给出一个试验方法，如图 1 所示。
|
||||
|
||||
<center>
|
||||
<img src="./img/CircleSquare.png">
|
||||
|
||||
图 1 正方形与内切圆
|
||||
</center>
|
||||
|
||||
中学几何的知识告诉我们：
|
||||
- 圆的面积 $S_c=\pi r^2$；
|
||||
- 正方形的面积 $S_s=2r \times 2r=4r^2$；
|
||||
- 那么圆的面积除以正方形的面积 $S_c/S_s=\pi r^2 / 4r^2=\pi/4$；
|
||||
- 所以 $\pi = 4S_c/S_s$。
|
||||
|
||||
如果是古人的话，可能会在图 1 上平铺一层细沙，然后称出在圆内的细沙的重量，再称出方形内细沙的重量，以此来模拟面积计算。现在我们借助计算机的力量，在图 1 内随机打点，然后统计出圆内的点的数量和方形内点的数量，也同样可以模拟出面积计算。
|
||||
|
||||
#### 实现代码
|
||||
|
||||
【代码位置：CirclePi_1.py】
|
||||
|
||||
```python
|
||||
# 随机投点
|
||||
def put_points(ax, num_total_points):
|
||||
ax.axis('equal')
|
||||
data = np.random.uniform(-1, 1, size=(num_total_points, 2)) # 在正方形内生成随机点
|
||||
r = np.sqrt(np.power(data[:,0], 2) + np.power(data[:,1], 2)) # 计算每个点到中心的距离
|
||||
num_in_circle = 0 # 统计在圆内的点数
|
||||
for i, point in enumerate(data): # 绘图
|
||||
if (r[i] < 1):
|
||||
num_in_circle += 1 # 计数
|
||||
ax.scatter(point[0], point[1], s=1, c='r')
|
||||
else:
|
||||
ax.scatter(point[0], point[1], s=1, c='b')
|
||||
# 计算 pi 值
|
||||
title = str.format("n={0},$\pi$={1}",num_total_points, num_in_circle/num_total_points*4)
|
||||
ax.set_title(title)
|
||||
draw_circle(ax, 1, 0, 0)
|
||||
ax.grid()
|
||||
```
|
||||
|
||||
为了避免随机性，可以使用相同的点数做多次试验并求平均，或者使用不同的点数做多次试验求平均。在此我们使用了第二种方法，绘制出图 2 所示的结果。
|
||||
|
||||
<center>
|
||||
<img src="./img/CirclePi-4.png">
|
||||
|
||||
图 2 随机投点估算圆周率
|
||||
</center>
|
||||
|
||||
在本次试验中，随机点数量与估算出的 $\pi$ 值对应关系如下：
|
||||
|
||||
- 100个点：3.04
|
||||
- 200个点：3.04
|
||||
- 500个点：3.136
|
||||
- 1000个点：3.18
|
||||
- 平均：3.099
|
||||
|
||||
人们通常以为点数越多，估算出的值越准确。从趋势上来看确实如此，但也不能避免单次试验的结果都带来的偏差。因此，在图 3 的试验中，我们分别使用 1000 到 20,000 个点来做估算，并且每个点数上都重复了 100 次取平均，最后得到了误差的变化图。
|
||||
|
||||
【代码位置：CirclePi_2.py】
|
||||
|
||||
<center>
|
||||
<img src="./img/CirclePi-error.png" width="500">
|
||||
|
||||
图 3 随机点计算圆周率的误差变化
|
||||
</center>
|
||||
|
||||
从横坐标上看，投点数量越多，偏差和方差都会随之减小。从最终的平均结果来看，3.1418 已经可以代表这种方法的精度了。当然，读者还可以用更多的点数来验证，不过单次试验的结果总会令人失望的。
|
||||
|
||||
对于估算不规则图形面积的问题，用这种方法同样可以解决。
|
||||
|
||||
### 2 提出问题二
|
||||
|
||||
#### 计算定积分
|
||||
|
||||
对于复杂函数，如何计算其定积分？
|
||||
|
||||
在图 4 中，显示了三个函数的图形，从左到右分别是：
|
||||
|
||||
- $f_1(x) = \frac{1}{x^2}$
|
||||
- $f_2(x) = \sin(x)$
|
||||
- $f_3(x) = 0.4x^2 + 0.3x\sin(15x) + 0.01\cos(50x)-0.3$
|
||||
|
||||
前两个简单函数，可以得知：
|
||||
|
||||
$$
|
||||
\int_{0.2}^1 \frac{1}{x^2} dx = 4, \quad \int_0^{\pi} \sin(x) dx = 2
|
||||
$$
|
||||
|
||||
那么第三个函数就要稍微费劲儿些了。
|
||||
|
||||
<center>
|
||||
<img src="./img/Integral.png">
|
||||
|
||||
图 4 计算函数积分
|
||||
</center>
|
||||
|
||||
#### 用蒙特卡洛方法估算定积分值
|
||||
|
||||
从问题一中，我们学到了随机投点的方法，在此仍然可以使用，但是由于函数值域区间的变化，写出代码来有些复杂（但还是可以完成的）。所以这一次介绍一个不同的方法，见图 5。
|
||||
|
||||
<center>
|
||||
<img src="./img/Integral-2.png">
|
||||
|
||||
图 5 计算函数积分
|
||||
</center>
|
||||
|
||||
假设把定积分的上下限 [a,b] 分成 N 等份，形成 N 个黄色的矩形，则其宽度为 $\frac{b-a}{N}$，其高度取矩形宽度的中点位置的 $f(x)$ 值，那么每个矩形的面积 $S_{x_i}=\frac{b-a}{N}f(x_i)$，于是曲线下的面积（即积分）可以近似为：
|
||||
|
||||
$$
|
||||
S = \sum_{i=1}^N S_{x_i} = \frac{b-a}{N} \sum_{i=1}^N f(x_i) \tag{1}
|
||||
$$
|
||||
|
||||
如果 $N$ 足够大的话，$S$ 的值会非常接近于数学积分值。
|
||||
|
||||
其理论基础是，随机变量 $X$ 的数学期望为：
|
||||
|
||||
$$
|
||||
\mathbb E[X]=\int_{-\infty}^{\infty} pdf(x)·x \ dx \tag{2}
|
||||
$$
|
||||
|
||||
其中 $pdf(x)$ 是概率分布函数，变成离散型公式就是：
|
||||
|
||||
$$
|
||||
S = \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{pdf(x_i)} \tag{3}
|
||||
$$
|
||||
|
||||
如果我们使用均匀分布，则 $pdf(X_i)=\frac{1}{b-a}$，是一个常量，所以式 3 可以变形为式 1。
|
||||
|
||||
#### 实现代码
|
||||
|
||||
【代码位置：Integral.py】
|
||||
|
||||
```Python
|
||||
def f1(x): # 函数 1
|
||||
y = 1/(x*x)
|
||||
return y
|
||||
|
||||
def f2(x): # 函数 2
|
||||
y = np.sin(x)
|
||||
return y
|
||||
|
||||
def f3(x): # 函数 3
|
||||
y = 0.4 * x * x + 0.3 * x * np.sin(15*x) + 0.01 * np.cos(50*x) - 0.3
|
||||
return y
|
||||
|
||||
def integral(f, a, b, n): # 计算积分
|
||||
v = 0
|
||||
repeat = 10 # 重复 10 次取平均
|
||||
for i in range(repeat):
|
||||
x = np.random.uniform(a, b, size=(n, 1)) # 随机生成 x
|
||||
y = f(x)
|
||||
v += np.sum(y) / n * (b-a) # 按式 1 计算
|
||||
return v/repeat
|
||||
```
|
||||
最后得到的结果是：
|
||||
```
|
||||
S1 = 4.008150279723373
|
||||
S2 = 2.0017614436591833
|
||||
S3 = -0.15136664866974026
|
||||
```
|
||||
|
||||
### 3 蒙特卡洛方法
|
||||
|
||||
#### 什么是蒙特卡洛方法
|
||||
|
||||
蒙特卡罗方法也称统计模拟方法，是一类随机算法的统称，是在1940年代中期由于科学技术的发展和电子计算机的发明，而提出的一种以概率统计理论为指导的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。
|
||||
|
||||
20世纪40年代，在冯·诺伊曼，斯塔尼斯拉夫·乌拉姆和尼古拉斯·梅特罗波利斯在洛斯阿拉莫斯国家实验室为核武器计划工作时，发明了蒙特卡罗方法。因为乌拉姆的叔叔经常在摩纳哥的蒙特卡洛赌场输钱得名，而蒙特卡罗方法正是以概率为基础的方法。与它对应的是确定性算法。
|
||||
|
||||
蒙特卡洛方法是一种近似推断的方法，通过采样大量样本的方法来求解期望、均值、面积、积分等问题。蒙特卡洛对某一种分布的采样方法有直接采样、接受拒绝采样与重要性采样三种，直接采样最简单，但是需要已知累积分布的形式。接受拒绝采样与重要性采样适用于原分布未知的情况，这两种方法都是给出一个提议分布，不同的是接受拒绝采样对不满足原分布的粒子予以拒绝，而重要性采样则是给予每个粒子不同的权重，大家可以根据不同的场景使用这三种方法中的一种进行采样。
|
||||
|
||||
读者不要单纯地认为点越多越精确，这是有一个前提的。总体而言，不精确是因为我们的点还没有到某个量级。一般情况下，蒙特卡罗算法的特点是，采样越多，越近似最优解，而永远不是最优解。
|
||||
|
||||
其优点比较明显，尤其对于具有统计性质的问题可以直接进行解决，对于连续性的问题也不必进行离散化处理。其缺点也很显然，对于确定性问题转化成随机性问题做的估值处理，丧失精确性，得到一个接近准确的N值也不太容易。
|
||||
|
||||
|
||||
#### 直接采样
|
||||
|
||||
#### 接受拒绝采样
|
||||
|
||||
#### 重要性采样
|
||||
|
||||
#### 蒙特卡洛方法分类
|
||||
|
||||
通常蒙特卡洛方法可以粗略地分成两类：
|
||||
|
||||
- 所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接**模拟随机的过程**
|
||||
例如在核物理研究中，分析中子在反应堆中的传输过程。中子与原子核作用受到量子力学规律的制约，人们只能知道它们相互作用发生的概率，却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向。
|
||||
科学家依据其概率进行随机抽样得到裂变位置、速度和方向，这样模拟大量中子的行为后，经过统计就能获得中子传输的范围，作为反应堆设计的依据。
|
||||
|
||||
- 所求解问题可以转化为某种**随机分布的特征数**
|
||||
比如随机事件出现的概率，或者随机变量的期望值。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解这种方法多用于求解复杂的多维积分问题。
|
||||
我们在上面举的例子是一个一维的积分问题，积分结果是曲线下的面积；当然可以再增加一维随机采样，变成二维积分问题，结果为曲面下的体积，比一维积分更有实际作用。
|
||||
|
||||
|
||||
### 参考资料
|
||||
|
||||
https://new.qq.com/omn/20220314/20220314A09IY600.html
|
||||
@@ -0,0 +1,63 @@
|
||||
|
||||
|
||||
### 7 如何确定采样的次数
|
||||
|
||||
<center>
|
||||
<img src="./img/Sampling-RMSE.png" width="500">
|
||||
|
||||
图 2
|
||||
</center>
|
||||
|
||||
|
||||
|
||||
具体地说就是如何确定算法中的幕数 Episodes 的数值。
|
||||
|
||||
根据问题的复杂程度不同，幕数必然会不同。但到目前为止，没有人从理论层面研究过这个问题，所以有一些偏实践的方法，供大家参考。
|
||||
|
||||
#### 试探
|
||||
|
||||
先用比较小的数值，比如 100，去做几次尝试，如果发现几次尝试的结果之间有很大的方差，就增加到 1000 再试试 ...... 以此类推，也许到 10 万时才能相对稳定。
|
||||
|
||||
方差公式为：
|
||||
|
||||
$$
|
||||
\sigma^2=\frac{1}{n}\sum_{i=1}^n (V_i-\mu)^2 \tag{2}
|
||||
$$
|
||||
|
||||
比如，我们只关注 Rest 状态的价值函数值，运行 3 次得到 $V_1,V_2,V_3$ 的值，而 $\mu=\frac{1}{3}(V_1+V_2+V_3)$ 是均值：
|
||||
|
||||
1. 当 episodes=100 时，运行 3 次，结果分别是 [1.2, 2.8, 2.0]，相差非常大，$\sigma^2=0.427$；
|
||||
2. 设置 episodes=1000，运行 3 次，结果分别是 [0.5, 1.2, 0.9]，方差减小了，$\sigma^2=0.082$，但还不够好；
|
||||
3. 设置 episodes=10000，运行 3 次，结果分别是 [0.86, 0.75, 0.82]，$\sigma^2=0.002$，如果这个方差到了你的心理预期，就可以结束了，否则可以再增加 episodes 的次数。
|
||||
|
||||
#### 比较
|
||||
|
||||
先看一个增量计算平均值的公式：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
V_{n+1} &= \frac{1}{n+1} \sum_{i=1}^{n+1} G_i=\frac{1}{n+1}(\sum_{i=1}^{n} G_i+G_{n+1})
|
||||
\\\\
|
||||
&= \frac{1}{n+1}(G_{n+1}+nV_n)= \frac{1}{n+1}(G_{n+1}+n V_n+ V_n-V_n)
|
||||
\\\\
|
||||
&=V_n + \frac{1}{n+1}(G_{n+1}-V_n)
|
||||
\end{aligned}
|
||||
\tag{3}
|
||||
$$
|
||||
|
||||
式 3 表达的意思是，n+1 幕时 $G_{n+1}$ 的期望值 $V_{n+1}$，等于 $n$ 幕时的 $G_n$ 期望值 $V_{n}$，再加上 n+1 幕时的 $G_{n+1}$ 与 $V_n$ 的差值除以 (n+1)。
|
||||
|
||||
把式 3 变形得到：
|
||||
|
||||
$$
|
||||
V_{n+1} - V_{n} = \frac{1}{n+1}(G_{n+1}-V_n)
|
||||
\tag{4}
|
||||
$$
|
||||
|
||||
|
||||
在第 n+1 幕时，先计算出式 4 的等号后面的部分，检查这个值是否足够小（比如小于 1e-2），就可以认为已经收敛了。当然，这里的步长可以不是 1 幕，而是 10 幕或者 100 幕才会做一次检查。
|
||||
|
||||
### 问题与讨论
|
||||
|
||||
1. 请使用试探法来找到比较理想的分幕次数。
|
||||
2. 请使用比较法来找到比较理想的分幕次数。
|
||||
@@ -0,0 +1,145 @@
|
||||
|
||||
## 每次访问型 MC
|
||||
|
||||
|
||||
|
||||
|
||||
### 算法改进
|
||||
|
||||
在上一节中，我们借助回报 $G$ 的定义
|
||||
|
||||
$$
|
||||
G_t = R_{t+1}+\gamma R_{t+2}+\gamma^2 R_{t+3}+ \cdots +\gamma^{T-t-1} R_{T} \tag{1}
|
||||
$$
|
||||
|
||||
以及价值函数 $V$ 的定义
|
||||
|
||||
$$
|
||||
V_t(s) = \mathbb E [G_t | S_t = s]
|
||||
\tag{2}
|
||||
$$
|
||||
|
||||
初步计算出了安全驾驶问题的各个状态的价值函数。回忆其过程如下：
|
||||
|
||||
1. 使用蒙特卡洛采样，每次采样都需要指定一个初始状态，然后在幕内循环，直到终止状态。
|
||||
2. 然后根据式 1 开始“从头”计算这个初始状态的回报 $G$ 值。
|
||||
3. 进行下一次采样，再计算 $G$ 值。
|
||||
4. 最后求平均（数学期望）得到 $V$。
|
||||
|
||||
聪明的读者可能会发现一个问题：如果不“从头”开始，而是从第二个、第三个状态开始计算，是不是就能在一次采样中就可以得到很多状态的 G 值呢？
|
||||
|
||||
<center>
|
||||
<img src="./img/MC-1.png" width="600">
|
||||
|
||||
图 1
|
||||
</center>
|
||||
|
||||
如图 1 所示，从 $S_1$ 开始一幕的采样，到 $S_T$ 为止结束，得到 $R$ 的序列后：
|
||||
- 固然可以从 $R_2$ 开始计算出 $G_1$。
|
||||
- 但是如果从 $R_3$ 开始，不就能计算出 $G_2$ 了吗？
|
||||
- 同理，还可以计算出这一采样序列中的任意的 $G_t$ 出来。
|
||||
|
||||
这样，利用一次采样结果可以计算出很多状态的 $G$ 值，会大幅提高算法的效率。
|
||||
|
||||
|
||||
### 算法描述
|
||||
|
||||
下面的伪代码中，$\leftarrow$ 表示赋值，$\Leftarrow$ 表示追加列表。
|
||||
|
||||
---
|
||||
|
||||
输入：起始状态$S,Episodes,\gamma$
|
||||
初始化：$G_{value}[S] \leftarrow 0, G_{count}[S] \leftarrow 0$
|
||||
多幕 $Episodes$ 循环：
|
||||
列表 $T = [\ ] $ 用于存储序列数据 $(S,R)$
|
||||
获得状态 $S$ 的奖励值 $R$
|
||||
$T \Leftarrow (S,R)$
|
||||
幕内循环直到终止状态：
|
||||
从 $S$ 根据状态转移概率得到 $S',R'$ 以及终止标志
|
||||
$T \Leftarrow S',R'$
|
||||
$S \leftarrow S'$
|
||||
对 $T$ 从后向前遍历, $t=T-1,T-2,...,0$
|
||||
从 $T$ 中取出 $S_t,R_t$
|
||||
$G \leftarrow \gamma G+R_t$
|
||||
$G_{value}[S_t] \leftarrow G_{value}[S_t]+G$
|
||||
$G_{count}[S_t] \leftarrow G_{value}[S_t]+1$
|
||||
$V[S] \leftarrow G_{value}[S] / G_{count}[S]$
|
||||
输出：$V[S]$
|
||||
|
||||
---
|
||||
|
||||
### 算法说明
|
||||
|
||||
<center>
|
||||
<img src="./img/MC-2.png">
|
||||
|
||||
图 2
|
||||
</center>
|
||||
|
||||
|
||||
### 算法实现
|
||||
|
||||
```Python
|
||||
# 反向计算G值，记录每个状态的G值，每次访问型
|
||||
def MC_Sampling_Reverse(dataModel, start_state, episodes, gamma):
|
||||
V_value_count_pair = np.zeros((dataModel.num_states, 2)) # state[total value, count of g]
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = [] # 按顺序 t 保留一幕内的采样序列
|
||||
trajectory.append((start_state.value, dataModel.get_reward(start_state)))
|
||||
curr_s = start_state
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, r, is_end = dataModel.step(curr_s)
|
||||
trajectory.append((next_s.value, r))
|
||||
curr_s = next_s
|
||||
#endwhile
|
||||
G = 0
|
||||
# 从后向前遍历
|
||||
for t in range(len(trajectory)-1, -1, -1):
|
||||
s, r = trajectory[t]
|
||||
G = gamma * G + r
|
||||
V_value_count_pair[s, 0] += G # 累积总和
|
||||
V_value_count_pair[s, 1] += 1 # 累计次数
|
||||
#endfor
|
||||
#endfor
|
||||
V = V_value_count_pair[:,0] / V_value_count_pair[:,1] # 计算平均值
|
||||
return V
|
||||
```
|
||||
|
||||
|
||||
|
||||
<center>
|
||||
<img src="./img/MC-2-RMSE.png" width="500">
|
||||
|
||||
图 2
|
||||
</center>
|
||||
|
||||
### 运行结果
|
||||
|
||||
表 $\gamma=1$ 时的状态值计算结果比较
|
||||
|
||||
|状态|原始算法|改进算法|准确值|
|
||||
|-|-:|-:|-:|-:|
|
||||
|出发 Start| 1.09|1.07|1.03|
|
||||
|正常行驶 Normal| 1.64|1.77|1.72|
|
||||
|礼让行人 Pedestrians| 2.74|2.75|2.72|
|
||||
|闹市减速 DownSpeed| 2.99|3.18|3.02|
|
||||
|超速行驶 ExceedSpeed| -5.81|-5.06|-5.17|
|
||||
|路口闯灯 RedLight| -6.72|-6.72|-6.73|
|
||||
|小区减速 LowSpeed| 6.00|6.00|6.00|
|
||||
|拨打电话 MobilePhone| -2.32|-2.40|-2.40|
|
||||
|发生事故 Crash| -1.00|-1.00|-1.00|
|
||||
|安全抵达 Goal| +5.00|5.00|5.00|
|
||||
|终止 End| 0.00| 0.00|0.00|
|
||||
|**误差 RMSE**|**0.042**|**0.034**||
|
||||
|
||||
可以看到，改进的算法在速度上和精度上都比原始算法要好。最后一行不是状态值，是 RMSE 的误差值，原始算法误差为 0.042，改进算法为 0.034，越小越好。
|
||||
|
||||
从性能上看，原始算法对每个状态做了 10000 次采样，相当于一共 $11 \times 10000=110000$ 次采样。改进算法对所有状态（混合）一共做了 50000 次采样。
|
||||
|
||||
|
||||
|
||||
### 参考资料
|
||||
|
||||
https://new.qq.com/omn/20220314/20220314A09IY600.html
|
||||
@@ -0,0 +1,7 @@
|
||||
|
||||
### 批量更新算法
|
||||
|
||||
|
||||
$$
|
||||
V(s) = V(s) + \alpha (G_t - V)
|
||||
$$
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 72 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 53 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 9.4 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 8.0 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 49 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 38 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 22 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 59 KiB |
@@ -0,0 +1,45 @@
|
||||
import random
|
||||
from turtle import color
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
# 画圆
|
||||
def draw_circle(ax,r,x,y):
|
||||
# 点的横坐标为a
|
||||
a = np.arange(x-r,x+r,0.000001)
|
||||
# 点的纵坐标为b
|
||||
b = np.sqrt(np.power(r,2)-np.power((a-x),2))+y
|
||||
ax.plot(a,b,color='g',linestyle='-')
|
||||
ax.plot(a,-b,color='g',linestyle='-')
|
||||
|
||||
# 随机投点
|
||||
def put_points(ax, num_total_points):
|
||||
ax.axis('equal')
|
||||
data = np.random.uniform(-1, 1, size=(num_total_points, 2)) # 在正方形内生成随机点
|
||||
r = np.sqrt(np.power(data[:,0], 2) + np.power(data[:,1], 2)) # 计算每个点到中心的距离
|
||||
num_in_circle = 0 # 统计在圆内的点数
|
||||
for i, point in enumerate(data): # 绘图
|
||||
if (r[i] < 1):
|
||||
num_in_circle += 1 # 计数
|
||||
ax.plot(point[0], point[1], 'o', markersize=1, c='r')
|
||||
else:
|
||||
ax.plot(point[0], point[1], 'o', markersize=1, c='b')
|
||||
# 计算 pi 值
|
||||
title = str.format("n={0},$\pi$={1}",num_total_points, num_in_circle/num_total_points*4)
|
||||
ax.set_title(title)
|
||||
draw_circle(ax, 1, 0, 0)
|
||||
ax.grid()
|
||||
|
||||
if __name__=="__main__":
|
||||
fig = plt.figure()
|
||||
|
||||
ax = fig.add_subplot(141)
|
||||
put_points(ax, 100)
|
||||
ax = fig.add_subplot(142)
|
||||
put_points(ax, 200)
|
||||
ax = fig.add_subplot(143)
|
||||
put_points(ax, 500)
|
||||
ax = fig.add_subplot(144)
|
||||
put_points(ax, 1000)
|
||||
|
||||
plt.show()
|
||||
@@ -0,0 +1,31 @@
|
||||
import tqdm
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
|
||||
|
||||
# 随机铺点
|
||||
def put_points(num_total_points):
|
||||
data = np.random.uniform(-1, 1, size=(num_total_points, 2)) # 在正方形内生成随机点
|
||||
r = np.sqrt(np.power(data[:,0], 2) + np.power(data[:,1], 2)) # 计算每个点到中心的距离
|
||||
r[r<=1]=1
|
||||
r[r>1]=0
|
||||
num_in_circle = np.sum(r)
|
||||
pi = num_in_circle/num_total_points*4
|
||||
return pi
|
||||
|
||||
if __name__=="__main__":
|
||||
pis = []
|
||||
for n in tqdm.trange(1000,20000,100):
|
||||
pi = 0
|
||||
for j in range(100):
|
||||
pi += put_points(n)
|
||||
pis.append(pi/100)
|
||||
|
||||
plt.grid()
|
||||
plt.plot(pis)
|
||||
plt.plot([0,200],[3.14159265,3.14159265])
|
||||
plt.title(str.format("average={0}", np.mean(pis)))
|
||||
plt.show()
|
||||
|
||||
# print(put_points(1000000))
|
||||
@@ -0,0 +1,81 @@
|
||||
|
||||
from asyncio import create_subprocess_shell
|
||||
import numpy as np
|
||||
from enum import Enum
|
||||
|
||||
# 状态
|
||||
class States(Enum):
|
||||
Start = 0 # 出发
|
||||
Normal = 1 # 正常行驶
|
||||
Pedestrians = 2 # 礼让行人
|
||||
DownSpeed = 3 # 闹市减速
|
||||
ExceedSpeed = 4 # 超速行驶
|
||||
RedLight = 5 # 路口闯灯
|
||||
LowSpeed = 6 # 小区减速
|
||||
MobilePhone = 7 # 拨打手机
|
||||
Crash = 8 # 发生事故
|
||||
Goal = 9 # 安全抵达
|
||||
End = 10 # 结束
|
||||
|
||||
# 状态转移概率
|
||||
P = np.array(
|
||||
[
|
||||
[0.0, 0.9, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.2, 0.1, 0.1, 0.1, 0.3, 0.1, 0.0, 0.1, 0.0],
|
||||
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.7, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.2, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.5, 0.0, 0.0],
|
||||
[0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
|
||||
[0.0, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
|
||||
]
|
||||
)
|
||||
|
||||
# 奖励向量
|
||||
# |出发|正常行驶|礼让行人|闹市减速|超速行驶|路口闯灯|小区减速|拨打手机|发生事故|安全抵达|结束|
|
||||
R = [0, 0, +1, +1, -3, -6, +1, -3, -1, +5, 0]
|
||||
|
||||
class DataModel(object):
|
||||
def __init__(self):
|
||||
self.P = P # 状态转移矩阵
|
||||
self.R = R # 奖励
|
||||
self.S = States # 状态集
|
||||
self.num_states = len(self.S) # 状态数量
|
||||
self.end_states = [self.S.End] # 终止状态集
|
||||
self.V_ground_truth = Matrix(self, 1)
|
||||
|
||||
# 判断给定状态是否为终止状态
|
||||
def is_end(self, s):
|
||||
if (s in self.end_states):
|
||||
return True
|
||||
return False
|
||||
|
||||
# 获得即时奖励，保留此函数可以为将来更复杂的奖励函数做准备
|
||||
def get_reward(self, s):
|
||||
return self.R[s.value]
|
||||
|
||||
def random_select_state(self):
|
||||
s = np.random.choice(self.num_states)
|
||||
return States(s)
|
||||
|
||||
# 根据转移概率前进一步，返回（下一个状态、即时奖励、是否为终止）
|
||||
def step(self, curr_s):
|
||||
next_s = np.random.choice(self.S, p=self.P[curr_s.value])
|
||||
return next_s, self.get_reward(next_s), self.is_end(next_s)
|
||||
|
||||
def Matrix(dataModel, gamma):
|
||||
num_state = dataModel.P.shape[0]
|
||||
I = np.eye(dataModel.num_states) * (1+1e-7)
|
||||
#I = np.eye(dataModel.num_states)
|
||||
tmp1 = I - gamma * dataModel.P
|
||||
tmp2 = np.linalg.inv(tmp1)
|
||||
vs = np.dot(tmp2, dataModel.R)
|
||||
return vs
|
||||
|
||||
if __name__=="__main__":
|
||||
dataModel = DataModel()
|
||||
v = Matrix(dataModel, 1.0)
|
||||
print(np.around(v,2))
|
||||
@@ -0,0 +1,101 @@
|
||||
|
||||
import tqdm
|
||||
import numpy as np
|
||||
import DriveDataModel as data
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
# 反向计算G值，记录每个状态的G值，每次访问型
|
||||
def MC_Sampling_Reverse_Checkpoint(dataModel, start_state, episodes, gamma, checkpoint):
|
||||
V = []
|
||||
V_value_count_pair = np.zeros((dataModel.num_states, 2)) # state[total value, count of g]
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = [] # 按顺序 t 保留一幕内的采样序列
|
||||
trajectory.append((start_state.value, dataModel.get_reward(start_state)))
|
||||
curr_s = start_state
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, r, is_end = dataModel.step(curr_s)
|
||||
trajectory.append((next_s.value, r))
|
||||
curr_s = next_s
|
||||
#endwhile
|
||||
G = 0
|
||||
# 从后向前遍历
|
||||
for t in range(len(trajectory)-1, -1, -1):
|
||||
s, r = trajectory[t]
|
||||
G = gamma * G + r
|
||||
V_value_count_pair[s, 0] += G # 累积总和
|
||||
V_value_count_pair[s, 1] += 1 # 累计次数
|
||||
#endfor
|
||||
if (episode+1)%checkpoint == 0:
|
||||
V.append(V_value_count_pair[:,0] / V_value_count_pair[:,1]) # 计算平均值
|
||||
#endfor
|
||||
V.append(V_value_count_pair[:,0] / V_value_count_pair[:,1]) # 计算平均值
|
||||
return V
|
||||
|
||||
# MC1的改进，反向计算G值，记录每个状态的G值，每次访问型
|
||||
def MC_Sampling_Reverse_FirstVisit(dataModel, start_state, episodes, gamma):
|
||||
V_value_count_pair = np.zeros((dataModel.num_states, 2)) # state[total value, count of g]
|
||||
V_value_count_pair[:,1] = 1 # 避免被除数为0
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = [] # 一幕内的采样序列
|
||||
curr_s = start_state
|
||||
trajectory.append((curr_s.value, dataModel.get_reward(start_state)))
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, r, is_end = dataModel.step(curr_s)
|
||||
#endif
|
||||
trajectory.append((next_s.value, r))
|
||||
curr_s = next_s
|
||||
#endwhile
|
||||
num_step = len(trajectory)
|
||||
G = 0
|
||||
first_visit = set()
|
||||
# 从后向前遍历
|
||||
for t in range(num_step-1, -1, -1):
|
||||
s, r = trajectory[t]
|
||||
G = gamma * G + r
|
||||
if (s in first_visit):
|
||||
continue
|
||||
V_value_count_pair[s, 0] += G # total value
|
||||
V_value_count_pair[s, 1] += 1 # count
|
||||
first_visit.add(s)
|
||||
#endfor
|
||||
#endfor
|
||||
V = V_value_count_pair[:,0] / V_value_count_pair[:,1]
|
||||
return V
|
||||
|
||||
def RMSE(a,b):
|
||||
err = np.sqrt(np.sum(np.square(a - b))/a.shape[0])
|
||||
return err
|
||||
|
||||
def test_once():
|
||||
episodes = 40000 # 计算 50000 次的试验的均值作为数学期望值
|
||||
gamma = 1 # 指定多个折扣因子做试验
|
||||
dataModel = data.DataModel()
|
||||
checkpoint = 100
|
||||
V = MC_Sampling_Reverse_Checkpoint(dataModel, dataModel.S.Start, episodes, gamma, checkpoint)
|
||||
print("gamma =", gamma)
|
||||
for s in dataModel.S:
|
||||
print(str.format("{0}:\t{1:.2f}", s.name, V[-1][s.value]))
|
||||
|
||||
errors = []
|
||||
for i in range(len(V)):
|
||||
err = RMSE(V[i], dataModel.V_ground_truth)
|
||||
errors.append(err)
|
||||
|
||||
return errors
|
||||
|
||||
if __name__=="__main__":
|
||||
ERRORS = []
|
||||
for i in range(10):
|
||||
errors = test_once()
|
||||
ERRORS.append(errors)
|
||||
|
||||
avg_E = np.mean(ERRORS, axis=0)
|
||||
plt.title(str.format("min RMSE={0}", np.min(avg_E)))
|
||||
plt.plot(avg_E)
|
||||
plt.grid()
|
||||
plt.show()
|
||||
@@ -0,0 +1,136 @@
|
||||
import multiprocessing as mp
|
||||
import tqdm
|
||||
import numpy as np
|
||||
import DriveDataModel as data
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# 多状态同时更新的蒙特卡洛采样
|
||||
# 注意输入V有初始状态
|
||||
# constant-alpha
|
||||
def MC_Batch_Update(dataModel, start_state, episodes, alpha, gamma, batch_size):
|
||||
V = np.zeros((dataModel.num_states))
|
||||
G_value_count_pair = np.zeros((dataModel.num_states,2)) # state[total value, count of g]
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = []
|
||||
if (start_state is None):
|
||||
curr_s = dataModel.random_select_state()
|
||||
else:
|
||||
curr_s = start_state
|
||||
|
||||
trajectory.append((curr_s.value, dataModel.get_reward(curr_s)))
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, R, is_end = dataModel.step(curr_s)
|
||||
#endif
|
||||
trajectory.append((next_s.value, R))
|
||||
curr_s = next_s
|
||||
|
||||
# calculate G_t
|
||||
num_step = len(trajectory)
|
||||
G = 0
|
||||
# 从后向前遍历
|
||||
for t in range(num_step-1, -1, -1):
|
||||
S, R = trajectory[t]
|
||||
G = gamma * G + R
|
||||
G_value_count_pair[S, 0] += G # total value
|
||||
G_value_count_pair[S, 1] += 1 # count
|
||||
|
||||
if ((episode+1)%batch_size == 0):
|
||||
G_batch = G_value_count_pair[:,0]/G_value_count_pair[:,1]
|
||||
V = V + alpha * (G_batch - V)
|
||||
G_value_count_pair[:,:] = 0
|
||||
#endfor
|
||||
#endfor
|
||||
return V
|
||||
|
||||
def MultiProcess(alphas):
|
||||
pool = mp.Pool(processes=4) # 指定合适的进程数量
|
||||
ERRORS = []
|
||||
results = []
|
||||
for alpha in alphas: # 遍历状态集中的每个状态作为起始状态
|
||||
results.append(pool.apply_async(test_once, args=(alpha,)))
|
||||
pool.close()
|
||||
pool.join()
|
||||
for i in range(len(alphas)):
|
||||
avg_e = results[i].get()
|
||||
ERRORS.append(avg_e)
|
||||
|
||||
return ERRORS # 4 个 alpha 对应的平均 error
|
||||
|
||||
# 每隔100幕保存一个中间结果
|
||||
def MC_Incremental_Update_Checkpoint(dataModel, start_state, episodes, alpha, gamma, checkpoint):
|
||||
Vs = []
|
||||
V = np.zeros((dataModel.num_states))
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = []
|
||||
if (start_state is None):
|
||||
curr_s = dataModel.random_select_state()
|
||||
else:
|
||||
curr_s = start_state
|
||||
|
||||
trajectory.append((curr_s.value, dataModel.get_reward(curr_s)))
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, R, is_end = dataModel.step(curr_s)
|
||||
#endif
|
||||
trajectory.append((next_s.value, R))
|
||||
curr_s = next_s
|
||||
|
||||
# calculate G_t
|
||||
num_step = len(trajectory)
|
||||
G = 0
|
||||
# 从后向前遍历
|
||||
for t in range(num_step-1, -1, -1):
|
||||
S, R = trajectory[t]
|
||||
G = gamma * G + R
|
||||
V[S] = V[S] + alpha * (G - V[S])
|
||||
#endfor
|
||||
if (episode+1)%checkpoint == 0:
|
||||
Vs.append(V.copy())
|
||||
#endfor
|
||||
return Vs
|
||||
|
||||
def RMSE(a,b):
|
||||
err = np.sqrt(np.sum(np.square(a - b))/a.shape[0])
|
||||
return err
|
||||
|
||||
# 针对一个alpha做10次的平均error
|
||||
def test_once(alpha):
|
||||
ERRORS = []
|
||||
for i in range(10):
|
||||
episodes = 20000 # 计算 10000 次的试验的均值作为数学期望值
|
||||
checkpoint = 100
|
||||
gamma = 1 # 指定多个折扣因子做试验
|
||||
dataModel = data.DataModel()
|
||||
errors = []
|
||||
Vs = MC_Incremental_Update_Checkpoint(dataModel, dataModel.S.Start, episodes, alpha, gamma, checkpoint)
|
||||
for result in Vs:
|
||||
err = RMSE(result, dataModel.V_ground_truth)
|
||||
errors.append(err)
|
||||
ERRORS.append(errors)
|
||||
avg_E = np.mean(ERRORS, axis=0)
|
||||
return avg_E
|
||||
|
||||
def test():
|
||||
alphas = [0.001,0.002,0.003,0.004]
|
||||
avgErros = MultiProcess(alphas)
|
||||
for i, avg_E in enumerate(avgErros):
|
||||
plt.plot(avg_E, label=str.format("alpha={0}, min RMSE={1:.3f}", alphas[i], np.min(avg_E)))
|
||||
plt.legend()
|
||||
plt.grid()
|
||||
plt.show()
|
||||
|
||||
def run():
|
||||
dataModel = data.DataModel()
|
||||
alpha = 0.02
|
||||
episodes = 20000
|
||||
gamma = 1
|
||||
batch_size = 100
|
||||
Vs = MC_Batch_Update(dataModel, dataModel.S.Start, episodes, alpha, gamma, batch_size)
|
||||
print(Vs)
|
||||
print(RMSE(Vs, dataModel.V_ground_truth))
|
||||
|
||||
if __name__=="__main__":
|
||||
run()
|
||||
@@ -0,0 +1,116 @@
|
||||
import multiprocessing as mp
|
||||
import tqdm
|
||||
import numpy as np
|
||||
import DriveDataModel as data
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# 多状态同时更新的蒙特卡洛采样
|
||||
# 注意输入V有初始状态
|
||||
# constant-alpha
|
||||
def MC_Incremental_Update(dataModel, start_state, episodes, alpha, gamma):
|
||||
V = np.zeros((dataModel.num_states))
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = []
|
||||
if (start_state is None):
|
||||
curr_s = dataModel.random_select_state()
|
||||
else:
|
||||
curr_s = start_state
|
||||
|
||||
trajectory.append((curr_s.value, dataModel.get_reward(curr_s)))
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, R, is_end = dataModel.step(curr_s)
|
||||
#endif
|
||||
trajectory.append((next_s.value, R))
|
||||
curr_s = next_s
|
||||
|
||||
# calculate G_t
|
||||
num_step = len(trajectory)
|
||||
G = 0
|
||||
# 从后向前遍历
|
||||
for t in range(num_step-1, -1, -1):
|
||||
S, R = trajectory[t]
|
||||
G = gamma * G + R
|
||||
V[S] = V[S] + alpha * (G - V[S])
|
||||
#endfor
|
||||
#endfor
|
||||
return V
|
||||
|
||||
def MultiProcess(alphas):
|
||||
pool = mp.Pool(processes=4) # 指定合适的进程数量
|
||||
ERRORS = []
|
||||
results = []
|
||||
for alpha in alphas: # 遍历状态集中的每个状态作为起始状态
|
||||
results.append(pool.apply_async(test_once, args=(alpha,)))
|
||||
pool.close()
|
||||
pool.join()
|
||||
for i in range(len(alphas)):
|
||||
avg_e = results[i].get()
|
||||
ERRORS.append(avg_e)
|
||||
|
||||
return ERRORS # 4 个 alpha 对应的平均 error
|
||||
|
||||
# 每隔100幕保存一个中间结果
|
||||
def MC_Incremental_Update_Checkpoint(dataModel, start_state, episodes, alpha, gamma, checkpoint):
|
||||
Vs = []
|
||||
V = np.zeros((dataModel.num_states))
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = []
|
||||
if (start_state is None):
|
||||
curr_s = dataModel.random_select_state()
|
||||
else:
|
||||
curr_s = start_state
|
||||
|
||||
trajectory.append((curr_s.value, dataModel.get_reward(curr_s)))
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, R, is_end = dataModel.step(curr_s)
|
||||
#endif
|
||||
trajectory.append((next_s.value, R))
|
||||
curr_s = next_s
|
||||
|
||||
# calculate G_t
|
||||
num_step = len(trajectory)
|
||||
G = 0
|
||||
# 从后向前遍历
|
||||
for t in range(num_step-1, -1, -1):
|
||||
S, R = trajectory[t]
|
||||
G = gamma * G + R
|
||||
V[S] = V[S] + alpha * (G - V[S])
|
||||
#endfor
|
||||
if (episode+1)%checkpoint == 0:
|
||||
Vs.append(V.copy())
|
||||
#endfor
|
||||
return Vs
|
||||
|
||||
def RMSE(a,b):
|
||||
err = np.sqrt(np.sum(np.square(a - b))/a.shape[0])
|
||||
return err
|
||||
|
||||
# 针对一个alpha做10次的平均error
|
||||
def test_once(alpha):
|
||||
ERRORS = []
|
||||
for i in range(2):
|
||||
episodes = 20000 # 计算 10000 次的试验的均值作为数学期望值
|
||||
checkpoint = 100
|
||||
gamma = 1 # 指定多个折扣因子做试验
|
||||
dataModel = data.DataModel()
|
||||
errors = []
|
||||
Vs = MC_Incremental_Update_Checkpoint(dataModel, dataModel.S.Start, episodes, alpha, gamma, checkpoint)
|
||||
for result in Vs:
|
||||
err = RMSE(result, dataModel.V_ground_truth)
|
||||
errors.append(err)
|
||||
ERRORS.append(errors)
|
||||
avg_E = np.mean(ERRORS, axis=0)
|
||||
return avg_E
|
||||
|
||||
if __name__=="__main__":
|
||||
alphas = [0.0005,0.001,0.002,0.005,0.01]
|
||||
avgErros = MultiProcess(alphas)
|
||||
for i, avg_E in enumerate(avgErros):
|
||||
plt.plot(avg_E, label=str.format("alpha={0}, min RMSE={1:.3f}", alphas[i], np.min(avg_E)))
|
||||
plt.legend()
|
||||
plt.grid()
|
||||
plt.show()
|
||||
@@ -0,0 +1,82 @@
|
||||
|
||||
import tqdm
|
||||
import numpy as np
|
||||
import DriveDataModel as data
|
||||
import time
|
||||
|
||||
# 反向计算G值，记录每个状态的G值，每次访问型
|
||||
def MC_Sampling_Reverse(dataModel, start_state, episodes, gamma):
|
||||
V_value_count_pair = np.zeros((dataModel.num_states, 2)) # state[total value, count of g]
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = [] # 按顺序 t 保留一幕内的采样序列
|
||||
trajectory.append((start_state.value, dataModel.get_reward(start_state)))
|
||||
curr_s = start_state
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, r, is_end = dataModel.step(curr_s)
|
||||
trajectory.append((next_s.value, r))
|
||||
curr_s = next_s
|
||||
#endwhile
|
||||
G = 0
|
||||
# 从后向前遍历
|
||||
for t in range(len(trajectory)-1, -1, -1):
|
||||
s, r = trajectory[t]
|
||||
G = gamma * G + r
|
||||
V_value_count_pair[s, 0] += G # 累积总和
|
||||
V_value_count_pair[s, 1] += 1 # 累计次数
|
||||
#endfor
|
||||
#endfor
|
||||
V = V_value_count_pair[:,0] / V_value_count_pair[:,1] # 计算平均值
|
||||
return V
|
||||
|
||||
# MC1的改进，反向计算G值，记录每个状态的G值，每次访问型
|
||||
def MC_Sampling_Reverse_FirstVisit(dataModel, start_state, episodes, gamma):
|
||||
V_value_count_pair = np.zeros((dataModel.num_states, 2)) # state[total value, count of g]
|
||||
V_value_count_pair[:,1] = 1 # 避免被除数为0
|
||||
for episode in tqdm.trange(episodes):
|
||||
trajectory = [] # 一幕内的采样序列
|
||||
curr_s = start_state
|
||||
trajectory.append((curr_s.value, dataModel.get_reward(start_state)))
|
||||
is_end = False
|
||||
while (is_end is False):
|
||||
# 从环境获得下一个状态和奖励
|
||||
next_s, r, is_end = dataModel.step(curr_s)
|
||||
#endif
|
||||
trajectory.append((next_s.value, r))
|
||||
curr_s = next_s
|
||||
#endwhile
|
||||
num_step = len(trajectory)
|
||||
G = 0
|
||||
first_visit = set()
|
||||
# 从后向前遍历
|
||||
for t in range(num_step-1, -1, -1):
|
||||
s, r = trajectory[t]
|
||||
G = gamma * G + r
|
||||
if (s in first_visit):
|
||||
continue
|
||||
V_value_count_pair[s, 0] += G # total value
|
||||
V_value_count_pair[s, 1] += 1 # count
|
||||
first_visit.add(s)
|
||||
#endfor
|
||||
#endfor
|
||||
V = V_value_count_pair[:,0] / V_value_count_pair[:,1]
|
||||
return V
|
||||
|
||||
def RMSE(a,b):
|
||||
err = np.sqrt(np.sum(np.square(a - b))/a.shape[0])
|
||||
return err
|
||||
|
||||
if __name__=="__main__":
|
||||
start = time.time()
|
||||
episodes = 50000 # 计算 10000 次的试验的均值作为数学期望值
|
||||
gammas = [0.5,0.9,1] # 指定多个折扣因子做试验
|
||||
dataModel = data.DataModel()
|
||||
for gamma in gammas:
|
||||
V = MC_Sampling_Reverse(dataModel, dataModel.S.Start, episodes, gamma)
|
||||
print("gamma =", gamma)
|
||||
for s in dataModel.S:
|
||||
print(str.format("{0}:\t{1:.2f}", s.name, V[s.value]))
|
||||
end = time.time()
|
||||
#print(end-start)
|
||||
#print(RMSE(V, dataModel.V_ground_truth))
|
||||
@@ -0,0 +1,59 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
def f1(x):
|
||||
y = 1/(x*x)
|
||||
return y
|
||||
|
||||
def f2(x):
|
||||
y = np.sin(x)
|
||||
return y
|
||||
|
||||
def f3(x):
|
||||
y = 0.4 * x * x + 0.3 * x * np.sin(15*x) + 0.01 * np.cos(50*x) - 0.3
|
||||
return y
|
||||
|
||||
def integral(f, a, b, n):
|
||||
v = 0
|
||||
repeat = 10
|
||||
for i in range(repeat):
|
||||
x = np.random.uniform(a, b, size=(n, 1))
|
||||
y = f(x)
|
||||
v += np.sum(y) / n * (b-a)
|
||||
return v/repeat
|
||||
|
||||
def show(ax, f, a, b, n, v):
|
||||
# 绘制函数曲线
|
||||
x = np.linspace(a, b, n)
|
||||
y = f(x)
|
||||
ax.set_title("integral="+str.format("{0:.2f}",v))
|
||||
ax.grid()
|
||||
ax.plot(x, y, c='g')
|
||||
|
||||
y_min = np.min(y)
|
||||
y_max = np.max(y)
|
||||
X = np.random.uniform(a, b, n)
|
||||
Y = np.random.uniform(y_min, y_max, n)
|
||||
for x,y in zip(X,Y):
|
||||
if y < f(x):
|
||||
ax.plot(x,y,'o',markersize=1,c='r')
|
||||
else:
|
||||
ax.plot(x,y,'o',markersize=1,c='b')
|
||||
|
||||
if __name__=="__main__":
|
||||
v1 = integral(f1, 0.2, 1, 10000)
|
||||
print("S1 =",v1)
|
||||
v2 = integral(f2, 0, 3.1416, 10000)
|
||||
print("S2 =",v2)
|
||||
v3 = integral(f3, 0, 1, 10000)
|
||||
print("S3 =",v3)
|
||||
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(131)
|
||||
show(ax, f1, 0.2, 1, 100, v1)
|
||||
ax = fig.add_subplot(132)
|
||||
show(ax, f2, 0, 3.1416, 100, v2)
|
||||
ax = fig.add_subplot(133)
|
||||
show(ax, f3, 0, 1, 100, v3)
|
||||
|
||||
plt.show()
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 24 KiB |
|
Before Width: | Height: | Size: 32 KiB After Width: | Height: | Size: 32 KiB |
|
Before Width: | Height: | Size: 34 KiB After Width: | Height: | Size: 34 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 38 KiB |
+29
-4
@@ -15,6 +15,7 @@ class States(Enum):
|
||||
# 奖励向量
|
||||
# [Game, Class1, Class2, Class3, Pass, Rest, End]
|
||||
Rewards = [-1, -2, -2, -2, 10, 1, 0]
|
||||
#Rewards = [-2, 0.5, 1, 1.5, 10, -2, 0]
|
||||
|
||||
# 状态转移概率
|
||||
P = np.array(
|
||||
@@ -25,20 +26,25 @@ P = np.array(
|
||||
[0.0, 0.0, 0.0, 0.0, 0.6, 0.4, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.2, 0.4, 0.4, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
|
||||
]
|
||||
)
|
||||
|
||||
ground_truth = [-22.54, -12.54, 1.46, 4.32, 10.00, 0.80, 0.00]
|
||||
ground_truth = [-22.543, -12.543, 1.457, 4.321, 10.00, 0.803, 0.00]
|
||||
a = [-22.5, -12.5, 1.5, 4.3, 10, 0.8, 0.00]
|
||||
|
||||
def RMSE(a,b):
|
||||
err = np.sqrt(np.sum(np.square(a - b))/a.shape[0])
|
||||
return err
|
||||
|
||||
class DataModel(object):
|
||||
def __init__(self):
|
||||
self.P = P # 状态转移矩阵
|
||||
self.R = Rewards # 奖励
|
||||
self.S = States # 状态集
|
||||
self.num_states = len(self.S) # 状态数量
|
||||
self.N = len(self.S) # 状态数量
|
||||
self.end_states = [self.S.End] # 终止状态集
|
||||
self.Y = ground_truth
|
||||
self.Y = SolveMatrix(self, 1)
|
||||
|
||||
# 判断给定状态是否为终止状态
|
||||
def is_end(self, s):
|
||||
@@ -54,3 +60,22 @@ class DataModel(object):
|
||||
def step(self, curr_s):
|
||||
next_s = np.random.choice(self.S, p=self.P[curr_s.value])
|
||||
return next_s, self.get_reward(next_s), self.is_end(next_s)
|
||||
|
||||
def SolveMatrix(dataModel, gamma):
|
||||
# 在对角矩阵上增加一个微小的值来解决奇异矩阵不可求逆的问题
|
||||
#I = np.eye(dataModel.N) * (1+1e-7)
|
||||
I = np.eye(dataModel.N)
|
||||
factor = I - gamma * dataModel.P
|
||||
inv_factor = np.linalg.inv(factor)
|
||||
vs = np.dot(inv_factor, dataModel.R)
|
||||
return vs
|
||||
|
||||
if __name__=="__main__":
|
||||
dataModel = DataModel()
|
||||
v = SolveMatrix(dataModel, 1.0)
|
||||
print(v)
|
||||
vv = np.around(v,3)
|
||||
for s in dataModel.S:
|
||||
print(str.format("{0}:\t{1}", s.name, vv[s.value]))
|
||||
print(RMSE(np.array(a), dataModel.Y))
|
||||
print(RMSE(np.array(ground_truth), dataModel.Y))
|
||||
@@ -0,0 +1,509 @@
|
||||
## 学生学习问题 - 贝尔曼方程
|
||||
|
||||
### 1 提出问题
|
||||
|
||||
在解决安全驾驶问题的过程中，我们学习了马尔可夫奖励过程、状态价值函数、蒙特卡洛采样法等强化学习的重要概念。蒙特卡洛法虽然是一种科学的方法，但是需要大量的采样才能得到**比较理想的结果**，并不能说是**准确的结果**。
|
||||
|
||||
蒙特卡洛法是针对**无模型**强化学习问题的，对于**有模型**的问题，我们有什么更好的方法可以解决吗？
|
||||
|
||||
在本节中，我们以在校大学生为例，描述学生学习（上课、复习、考试）的过程，来分析解决上述问题。
|
||||
|
||||
假设一门课只需要上三次课就可以结束，然后就可以通过考试而结课拿学分。当然，这其中也不是那么顺利的，学生可能会遇到各种挑战：
|
||||
|
||||
- 上课不专心听讲而去打手机游戏；
|
||||
- 学到一半的时候觉得这门课索然无味，中途退课；
|
||||
- 觉得离考试还远，不着急复习巩固知识，而是去休息；
|
||||
......
|
||||
|
||||
|
||||
### 2 建立模型
|
||||
|
||||
|
||||
|
||||
图 1 是一个有关学生的学习、考试等一些列状态的马尔可夫链，也可以叫做状态转移图，并且给每个状态都赋予了一个即时奖励值。
|
||||
|
||||
<center>
|
||||
<img src="./img/Student-2.png" width="500">
|
||||
|
||||
图 1 学习问题的状态转移概率图
|
||||
</center>
|
||||
|
||||
- Class 1,2,3
|
||||
上课/学习/复习状态，假设一门课是需要三次课的学习就可以结束。
|
||||
- 在上课 1 中，有 0.5 的概率跑到 Game 状态，在课堂上用手机偷偷摸摸打游戏，另外 0.5 的概率到上课 2 状态。
|
||||
- 在上课 2 中，有 0.2 的概率直接退课，另外 0.8 的概率到上课 3 状态。
|
||||
- 在上课 3 中，有 0.6 的概率去考试，另外 0.4 的概率以为考试还远，不着急准备呢，就去休息了。
|
||||
- Pass
|
||||
考试通过状态。考试结束后，以 100% 概率到结课状态。
|
||||
- Rest
|
||||
休息状态。休息结束后，发现学的内容都忘得差不多了，遂分别以不同的概率回到三次课的上课/学习/复习状态。
|
||||
在这里我们不区分上课和复习，如果是第一次到达 C1 状态，就认为是上课，否则就认为是复习。
|
||||
- Game
|
||||
娱乐/打游戏状态。在游戏状态中很大可能不能自拔，以 0.9 概率继续打游戏，只有 0.1 的概率幡然悔悟回到学习状态。
|
||||
- End
|
||||
结课/退课状态。进入此状态后将不再进行转移，或者是说以 100% 的概率转移到自己，叫做结束状态或者吸收状态。
|
||||
|
||||
表 1 中列出了状态转移矩阵，与租车问题中的矩阵形式相同。
|
||||
|
||||
表 1 状态转移矩阵
|
||||
|
||||
|P: 从$\rightarrow$到|Game|Class1|Class2|Class3|Pass|Rest|End|
|
||||
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|
||||
|**Game**|0.9|0.1||||||
|
||||
|**Class1**|0.5||0.5|||||
|
||||
|**Class2**||||0.8|||0.2|
|
||||
|**Class3**|||||0.6|0.4||
|
||||
|**Pass**|||||||1.0|
|
||||
|**Rest**||0.2|0.4|0.4||||
|
||||
|**End**|||||||1.0|
|
||||
|
||||
有的读者可能有个疑问：打游戏上瘾，从 Game 到 Game 有 0.9 的高概率，那么当该学生打游戏 2 小时后良心发现，转到学习状态的概率会不会大于 0.1 呢？
|
||||
|
||||
这是一个简单的平稳环境的马尔科夫链，如果考虑更复杂的情况，可以在 Game 状态下增加一个计数器：
|
||||
- 如果打游戏超过 1 小时了，则有 0.5 的概率回到学习状态；
|
||||
- 如果没超过 1 小时，则有 0.1 的概率回到学习状态。
|
||||
|
||||
|
||||
在学生学习的模型中有很多状态，如何确定某个状态比另一个状态好呢？或者说如何比较两个状态的好坏呢？因为从直觉上讲，学生在 C1,C2,C3 的状态明显要比 Game 状态好，但是如何能用数值的大小来体现这种好坏关系呢？
|
||||
|
||||
|
||||
在上一节中，已经有了分幕、奖励、回报的概念，这一节中，将会利用这些基础概念来定义每个**状态价值函数**，从而可以比较状态之间的好坏。
|
||||
|
||||
需要再次说明的是，在图 1 中，我们使用了**注重结果**的奖励定义方式，直接给每个状态赋值一个奖励，意味只要达到这个状态，就可以立刻获得标注出的奖励值，而不管是从哪条路径达到的。
|
||||
|
||||
|
||||
### 马尔可夫奖励过程中的贝尔曼方程
|
||||
|
||||
<center>
|
||||
<img src="./img/Bellman.png">
|
||||
|
||||
图 1 贝尔曼公式推导
|
||||
</center>
|
||||
|
||||
- 左图：在 $s_a$ 状态得到 $R(s)$ 的表达式。
|
||||
在使用**注重过程**的奖励函数定义方式时，从$s_a$ 转移到 $s_b,s_c$ 的过程中，分别可以得到 $r_1,r_2$ 的奖励，则 $s_a$ 的奖励函数定义为一种期望：$R(s_a)=\mathbb E[R_{Sa}|S_t=s_a] = p_1 \cdot r_1+p_2 \cdot r_2$。
|
||||
- 右图：在 $s_a$ 状态得到 $G_{t+1}$ 的表达式。
|
||||
当 $S_t=s$ 时，即在 $s_a$状态下，只能确定 $G_{t}=G_a$，不能确定$G_{t+1}$，因为不知道下一步会转移到哪个状态，是 $s_b$ 还是 $s_c$？
|
||||
所以，在 $s_a$ 状态时，$G_{t+1}$ 只能用转移概率（即 $p_1,p_2$）与下层状态 $s_b,s_c$ 的 $G$ 值（即 $G_a,G_b$）的乘积来表示，相当于在 $S_t=s_a$ 时，对 $G_{t+1}$ 求一次期望（带权重的平均值）：$G_{t+1}=\mathbb E[G_{b,c}|S_t=s_a]=(p_1 \cdot G_b|S_{t+1}=s_b)+(p_2 \cdot G_c|S_{t+1}=s_c)$。一旦确定到达 $s_b,s_c$ 状态后，$S_t=s_a$ 的条件就可以去掉了，分别用 $S_{t+1}=s_b,S_{t+1}=s_c$ 代替。
|
||||
|
||||
|
||||
做实例化推导之前，针对图 1，先给出一些必要的定义。
|
||||
|
||||
状态集定义：
|
||||
|
||||
$$
|
||||
s_a \in s, \ (s_b,s_c) \in s'
|
||||
$$
|
||||
|
||||
其中：$s$ 等同于 $S_t$，$s'$ 等同于 $S_{t+1}$。
|
||||
|
||||
由价值函数的定义：
|
||||
$$
|
||||
V(s)=\mathbb E [G_t|S_t=s] \tag{1}
|
||||
$$
|
||||
可以得到图 1 中状态 $S_a, S_b, S_c$ 的价值函数的实例化表示：
|
||||
$$
|
||||
\begin{aligned}
|
||||
V(s_a)&=\mathbb E [G_a|S_t=s_a] & (2.1)
|
||||
\\
|
||||
V(s_b)&=\mathbb E [G_b|S_{t+1}=s_b] & (2.2)
|
||||
\\
|
||||
V(s_c)&=\mathbb E [G_c|S_{t+1}=s_c] & (2.3)
|
||||
\end{aligned}
|
||||
\tag{2}
|
||||
$$
|
||||
|
||||
在本例中，如果 $V(s)=V(s_a)$，则 $V(s_b),V(s_c) \in V(s')$。
|
||||
|
||||
图 1 中状态转移概率的实例化表示：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
p_1 &= p(s_b|s_a)=P(s'|s), \ (s=s_a,s'=s_b) &(3.1)
|
||||
\\
|
||||
p_2 &= p(s_c|s_a)=P(s'|s), \ (s=s_a,s'=s_c) &(3.2)
|
||||
\end{aligned}
|
||||
\tag{3}
|
||||
$$
|
||||
|
||||
图 1 中奖励函数的实例化表示：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
r_1 &= r(s_a,s_b)=R(s,s'), \ (s=s_a,s'=s_b) &(4.1)
|
||||
\\
|
||||
r_2 &= r(s_a,s_c)=R(s,s'), \ (s=s_a,s'=s_c) &(4.2)
|
||||
\end{aligned}
|
||||
\tag{4}
|
||||
$$
|
||||
|
||||
该奖励函数的定义属于**注重过程**的定义方式，即定义在状态转移过程中。而此时状态 $S_a$ 的奖励为：
|
||||
$$
|
||||
R_{t+1}=p_1 \cdot r_1+p_2 \cdot r_2 \tag{5}
|
||||
$$
|
||||
|
||||
|
||||
推导
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
V(s)&=\mathbb E [G_t|S_t=s]
|
||||
\\
|
||||
&=\mathbb E[R_{t+1}+\gamma R_{t+2}+\gamma^2 R_{t+3}+\cdots|S_t=s]
|
||||
\\
|
||||
&=\mathbb E[R_{t+1}+\gamma (R_{t+2}+\gamma R_{t+3}+\cdots)|S_t=s]
|
||||
\\
|
||||
&=\mathbb E[R_{t+1}+\gamma G_{t+1}|S_t=s]
|
||||
\\
|
||||
&= \underbrace{ \mathbb E[R_{t+1}|S_t=s]}_A + \gamma \underbrace{\mathbb E[G_{t+1}|S_t=s]}_B
|
||||
\end{aligned}
|
||||
\tag{6}
|
||||
$$
|
||||
|
||||
式 6 的 $A$ 部分：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
A & = \mathbb E[R_{t+1}|S_t=s]
|
||||
\\
|
||||
({\footnotesize 实例化}\to) &= \mathbb E[R_{Sa}|S_t=s_a]
|
||||
\\
|
||||
({\footnotesize 图 1 左图} \to )&= p_1 \cdot r_1+p_2 \cdot r_2
|
||||
\\
|
||||
({\footnotesize 式3,4} \to) &=p(s_b|s_a) \cdot r(s_a,s_b)+p(s_c|s_a) \cdot r(s_a,s_c)
|
||||
\\
|
||||
({\footnotesize 抽象化} \to) &=\sum_{s'} P(s'|s) \cdot R(s,s') \to R(s)
|
||||
\end{aligned}
|
||||
\tag{7}
|
||||
$$
|
||||
|
||||
式 6 的 $B$ 部分：
|
||||
|
||||
首先要注意的一个问题是，B 不等于 $V(s')$，因为按式 6 的第一行，$V(s')=\mathbb E[G_{t+1}|S_{t+1}=s'] \ne \mathbb E[G_{t+1}|S_t=s]$。
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
B&=\mathbb E\big[G_{t+1}|S_t=s \big ]
|
||||
\\
|
||||
({\footnotesize 实例化} \to)&=\mathbb E \big[\mathbb E[G_{b,c}|S_t=s_a] \big ]
|
||||
\\
|
||||
({\footnotesize 图1右图} \to)&= \mathbb E\big[(p_1 \cdot G_{b}|S_{t+1}=s_b) + (p_2\cdot G_{c}|S_{t+1}=s_c)\big]
|
||||
\\
|
||||
({\footnotesize 期望加法变换} \to)&=\mathbb E\big[p_1\cdot G_{b}|S_{t+1}=s_b]+\mathbb E[p_2\cdot G_{c}|S_{t+1}=s_c\big]
|
||||
\\
|
||||
({\footnotesize 提出常数} p_1,p_2\to)&=p_1 \cdot \mathbb E[G_{b}|S_{t+1}=s_b]+ p_2 \cdot \mathbb E[G_{c}|S_{t+1}=s_c]
|
||||
\\
|
||||
({\footnotesize 式2} \to)&= p(s_b|s_a) \cdot V(s_b) + p(s_c|s_a) \cdot V(s_c)
|
||||
\\
|
||||
({\footnotesize 抽象化} \to)&= \sum_{s'} P(s'|s)V(s')
|
||||
\end{aligned}
|
||||
\tag{8}
|
||||
$$
|
||||
|
||||
所以式 6 最终为：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
V(s) &= \mathbb E[R_{t+1}|S_t=s] + \gamma \mathbb E[G_{t+1}|S_t=s]
|
||||
\\
|
||||
&=\sum_{s'} P(s'|s) R(s,s')+ \gamma \sum_{s'} P(s'|s)V(s') & (9.1)
|
||||
\\
|
||||
&=\sum_{s'} P(s'|s)[R(s,s')+\gamma V(s')] &(9.2)
|
||||
\\
|
||||
&= R(s)+ \gamma \sum_{s'} P(s'|s)V(s')=R_s+ \gamma \sum_{s'} P_{ss'}V(s') &(9.3)
|
||||
\end{aligned}
|
||||
\tag{9}
|
||||
$$
|
||||
|
||||
- 在**针对过程定义奖励函数**的问题中，使用式 9.2 比较方便，因为 $R(s,s')$ 是定义在从 $s\to s'$ 的转移过程上。这是 Richard S. Sutton and Andrew G. Barto 书中的写法。
|
||||
- 在**针对状态定义奖励函数**的问题中，使用式 9.3 比较方便，因为 $R(s)$ 是直接定义在状态 $s$ 上。这是 David Silver 课件中的写法。
|
||||
|
||||
|
||||
如果针对图 1，状态 $s_a$ 的价值函数实例计算公式为：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
V(s_a)&=(p_1 \cdot r_1 + p_2 \cdot r_2) + \gamma[p_1 \cdot V(s_b) + p_2 \cdot V(s_c)]
|
||||
\\
|
||||
&=R(s)+\gamma[p_1 \cdot V(s_b) + p_2 \cdot V(s_c)]
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
也就是说，一个状态 $s$ 的价值函数 $V(s)$ 由它的下游状态 $s'$ 的价值函数 $V(s')$ 和转移概率 $P(s,s')$ 以及转移过程中的奖励 $R(s,s')$ 构成。
|
||||
|
||||
Bellman Equation for MRP
|
||||
|
||||
<center>
|
||||
<img src="./img/student-3.png" width="500">
|
||||
|
||||
图 2
|
||||
</center>
|
||||
|
||||
图 2 中，每个状态下方都用括号表示了该状态的序号，比如 C1(1) 表示 $v_1$。以状态 C3 为例，根据式 9.3，可以得到其价值函数为：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
v_3&=R(C3)+\gamma[P_{C3,Pass} \cdot V(Pass) + P_{C3,Rest} \cdot V(Rest)]
|
||||
\\
|
||||
&=-2+ (0.6 v_4 + 0.4 v_5), &(\gamma=1)
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
同理可以得到其它所有状态的价值函数表达式，列出方程组如下：
|
||||
|
||||
$$
|
||||
\begin{cases}
|
||||
v_0=-1+0.9v_0+0.1v_1 & (10.1)
|
||||
\\
|
||||
v_1=-2+0.5v_0+0.5v_2 & (10.2)
|
||||
\\
|
||||
v_2=-2+0.8v_3+0.2v_6 & (10.3)
|
||||
\\
|
||||
v_3=-2+0.6v_4+0.4v_5 & (10.4)
|
||||
\\
|
||||
v_4=10+v_6 & (10.5)
|
||||
\\
|
||||
v_5=1+0.2v_1+0.4v_2+0.4v_3 & (10.6)
|
||||
\\
|
||||
v_6=0 & (10.7)
|
||||
\end{cases}
|
||||
\tag{10}
|
||||
$$
|
||||
|
||||
这是一个七元一次方程组，肯定有解。先简化式 10 中的各项，得到新的表达式：
|
||||
|
||||
$$
|
||||
\begin{cases}
|
||||
v_0=v_1-10 & (11.1)
|
||||
\\
|
||||
v_1=-2+0.5v_0+0.5v_2 & (11.2)
|
||||
\\
|
||||
v_2=0.8v_3-2 & (11.3)
|
||||
\\
|
||||
v_3=0.4v_5+4 & (11.4)
|
||||
\\
|
||||
v_4=10 & (11.5)
|
||||
\\
|
||||
v_5=1+0.2v_1+0.4v_2+0.4v_3 & (11.6)
|
||||
\\
|
||||
v_6=0 & (11.7)
|
||||
\end{cases}
|
||||
\tag{11}
|
||||
$$
|
||||
|
||||
将 $(11.1)(11.2)(11.3)(11.4)$ 都变成 $v_3$ 的表达式，带入$(11.6)$ 的两侧，可以得到：
|
||||
|
||||
$$
|
||||
2.5v_3-10=1+0.2(0.8v_3-16)+0.4(0.8v_3-2)+0.4v_3
|
||||
$$
|
||||
|
||||
得到：$v_3=4.321$
|
||||
|
||||
所以，最终的结果为：
|
||||
|
||||
$$
|
||||
\begin{cases}
|
||||
v_0=-22.543 \approx -22.5
|
||||
\\
|
||||
v_1=-12.543 \approx -12.5
|
||||
\\
|
||||
v_2=1.457 \approx 1.5
|
||||
\\
|
||||
v_3=4.321 \approx 4.3
|
||||
\\
|
||||
v_4=10
|
||||
\\
|
||||
v_5=0.803 \approx 0.8
|
||||
\\
|
||||
v_6=0
|
||||
\end{cases}
|
||||
\tag{12}
|
||||
$$
|
||||
|
||||
读者可以用式 12 的结果验证式 10 中的任意等式。
|
||||
|
||||
|
||||
### 矩阵法
|
||||
|
||||
观察式 10 方程组，可以把它变形为：
|
||||
|
||||
$$
|
||||
\begin{bmatrix}
|
||||
v_0
|
||||
\\
|
||||
v_1
|
||||
\\
|
||||
v_2
|
||||
\\
|
||||
v_3
|
||||
\\
|
||||
v_4
|
||||
\\
|
||||
v_5
|
||||
\\
|
||||
v_6
|
||||
\end{bmatrix}
|
||||
= \
|
||||
\begin{bmatrix}
|
||||
-1
|
||||
\\
|
||||
-2
|
||||
\\
|
||||
-2
|
||||
\\
|
||||
-2
|
||||
\\
|
||||
10
|
||||
\\
|
||||
1
|
||||
\\
|
||||
0
|
||||
\end{bmatrix}
|
||||
+\gamma
|
||||
\begin{bmatrix}
|
||||
0.9v_0+0.1v_1
|
||||
\\
|
||||
0.5v_0+0.5v_2
|
||||
\\
|
||||
0.8v_3+0.2v_6
|
||||
\\
|
||||
0.6v_4+0.4v_5
|
||||
\\
|
||||
v_6
|
||||
\\
|
||||
0.2v_1+0.4v_2+0.4v_3
|
||||
\\
|
||||
0
|
||||
\end{bmatrix}
|
||||
\tag{13}
|
||||
$$
|
||||
|
||||
|
||||
关于式 13：
|
||||
|
||||
- 等式左侧的部分，就是状态值的向量，可以写成 $V_s$。
|
||||
- 等式右侧的第一项，就是状态上的奖励值组成的向量，可以写成 $R_s$。
|
||||
- 等式右侧的第二个矩阵，又可以写成两个矩阵的乘积：
|
||||
|
||||
$$
|
||||
\begin{bmatrix}
|
||||
0.9 & 0.1 & 0 & 0 & 0 & 0 & 0
|
||||
\\
|
||||
0.5 & 0 & 0.5 & 0 & 0 & 0 & 0
|
||||
\\
|
||||
0 & 0 & 0 & 0.8 & 0 & 0 & 0.2
|
||||
\\
|
||||
0 & 0 & 0 & 0 & 0.6 & 0.4 & 0
|
||||
\\
|
||||
0 & 0 & 0 & 0 & 0 & 0 & 1.0
|
||||
\\
|
||||
0 & 0.2 & 0.4 & 0.4 & 0 & 0 & 0
|
||||
\\
|
||||
0 & 0 & 0 & 0 & 0 & 0 & 1.0
|
||||
\end{bmatrix}
|
||||
\begin{bmatrix}
|
||||
v_0
|
||||
\\
|
||||
v_1
|
||||
\\
|
||||
v_2
|
||||
\\
|
||||
v_3
|
||||
\\
|
||||
v_4
|
||||
\\
|
||||
v_5
|
||||
\\
|
||||
v_6
|
||||
\end{bmatrix}
|
||||
\tag{14}
|
||||
$$
|
||||
|
||||
其中，第一个矩阵就是该问题的状态转移矩阵 $P_{ss'}$，第二个矩阵是状态值向量 $V_s$，于是，式 9 可以变成：
|
||||
$$
|
||||
V_s = R_s+ \gamma P_{ss'}V_s \tag{15}
|
||||
$$
|
||||
|
||||
从式 9 直接看过来，式 15 等式右侧的 $V_s$ 应该是 $V_{s'}$ 才对，即 $V_s = R_s+\gamma P_{ss'}V_{s'}$。但是经过上述的实例化推导，读者可以理解所谓的 $V_{s'}$ 是在时间维度上的定义，表示下一步的状态；而在空间上，由于状态值一旦确定就不会变化，并没有 $s'$ 的概念。
|
||||
|
||||
比如：
|
||||
- 式 10.4，$v_3=-2+ (0.6 v_4 + 0.4 v_5)$ 中，$V_s=v_3,V_{s'}=\{v_4,v_5\}$，$v_5$ 是 $v_3$ 的后续状态。
|
||||
- 式 10.6，$v_5=1+0.2v_1+0.4v_2+0.4v_3$ 中，$V_s=v_5,V_{s'}=\{v_1,v_2,v_3\}$，$v_3$ 是 $v_5$ 的后续状态。
|
||||
|
||||
两者在不同的马尔可夫过程中互为后续状态，所以实际上并没有 $V_{s'}$ 的概念，或者说 $V_s$ 和 $V_{s'}$ 对于具体的状态实例有区别，对于状态向量组没有区别。
|
||||
|
||||
式 15 可以变形，并最终解出 $V_s$：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
V_s &= R_s+ \gamma P_{ss'}V_s
|
||||
\\
|
||||
V_s - \gamma P_{ss'}V_s &= R_s
|
||||
\\
|
||||
(I-\gamma P_{ss'})V_s&=R_s, &(I \ {\footnotesize 是对角矩阵})
|
||||
\\
|
||||
V_s&=(I-\gamma P_{ss'})^{-1}R_s
|
||||
\end{aligned}
|
||||
\tag{16}
|
||||
$$
|
||||
|
||||
式 16 中，等式右侧的值都是已知的，所以可以解出 $V_s$ 的数学解析解。
|
||||
|
||||
定义状态转移矩阵：
|
||||
|
||||
```python
|
||||
# 状态转移概率
|
||||
P = np.array(
|
||||
[ #Game Cl1 Cl2 Cl3 Pass Rest End
|
||||
[0.9, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.2],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.6, 0.4, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
|
||||
[0.0, 0.2, 0.4, 0.4, 0.0, 0.0, 0.0],
|
||||
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
|
||||
]
|
||||
)
|
||||
```
|
||||
|
||||
定义奖励函数值向量：
|
||||
```python
|
||||
# 奖励向量
|
||||
# [Game, Class1, Class2, Class3, Pass, Rest, End]
|
||||
Rewards = [-1, -2, -2, -2, 10, 1, 0]
|
||||
```
|
||||
代码
|
||||
|
||||
```python
|
||||
def SolveMatrix(dataModel, gamma):
|
||||
# 在对角矩阵上增加一个微小的值来解决奇异矩阵不可求逆的问题
|
||||
I = np.eye(dataModel.N) * (1+1e-7)
|
||||
# I = np.eye(dataModel.N) # 非奇异矩阵时使用此行代码以提高计算精度
|
||||
factor = I - gamma * dataModel.P
|
||||
inv_factor = np.linalg.inv(factor)
|
||||
vs = np.dot(inv_factor, dataModel.R)
|
||||
return vs
|
||||
```
|
||||
在定义状态转移矩阵时，右下角的值，即从 $S_{End} \to S_{End}$ 的转移概率，即可以写成 0.0，也可以写成 1.0，从强化学习的概念角度出发，都没有错。
|
||||
|
||||
与代码处理逻辑有关。
|
||||
|
||||
写成 0.0 时，np.random.choice(p=) 函数由于概率之和不为 1，所以函数调用出错。
|
||||
写成 1.0 时，由于矩阵的行列式为 0，是个奇异矩阵，不可求逆。此时可以在对角矩阵上增加一个微小的值来解决奇异矩阵不可求逆的问题，但是最终结果会有 1e-7 的误差，可以接受。
|
||||
|
||||
```
|
||||
[-22.54320988 -12.54320988 1.45679012 4.32098765 10. 0.80246914 0. ]
|
||||
Game: -22.543
|
||||
Class1: -12.543
|
||||
Class2: 1.457
|
||||
Class3: 4.321
|
||||
Pass: 10.0
|
||||
Rest: 0.802
|
||||
End: 0.0
|
||||
```
|
||||
|
||||
可能有读者好奇：用矩阵法得到的结果，与式 12 相比，哪一个更准确？
|
||||
|
||||
答案是：式 12 更准确。原因是用代码求矩阵的逆时，由于具体实现的问题有一些误差，否则的话两者应该完全相等。
|
||||
|
||||
### 迭代法（动态规划）
|
||||
@@ -8,3 +8,19 @@ draft里面将会以日记形式存放独立的小“故事”，每个“故事
|
||||
- img（存放图片）
|
||||
- src（存放代码）
|
||||
问题.md（正文内容）
|
||||
|
||||
|
||||
三门问题
|
||||
租车问题
|
||||
醉汉回家问题
|
||||
安全驾驶问题
|
||||
学生学习问题
|
||||
网格世界（A->A', B-B', 5x5)
|
||||
悬崖问题
|
||||
FrozenLake
|
||||
Car
|
||||
BlackJack
|
||||
肥皂泡
|
||||
有风的世界
|
||||
迷宫
|
||||
Gym里面的more
|
||||
|
||||
Reference in New Issue
Block a user