* update

* jjk

* update

* yi

* iuou

* update

* update

* update

* upodate

* Update 2-首次访问法与每次访问法.md

* update

* updae

* dfj

* update

* update

* update

* update

* up

* udate

* update

* update

* up

* update

* update

* update

* update

* update

* update

* update

* update

* update

* Update ValueIteration.png

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* update

* Update MC_111_ES.py

* update

* update

* update
This commit is contained in:
xiaowuhu 2022-08-02 12:04:35 +08:00 коммит произвёл GitHub
Родитель 96ed38c5ee
Коммит bd757f1a48
Не найден ключ, соответствующий данной подписи
Идентификатор ключа GPG: 4AEE18F83AFDEB23
136 изменённых файлов: 5318 добавлений и 681 удалений

Просмотреть файл

@ -0,0 +1,26 @@
|英文符号或缩写|英文单词或全拼|中文含义|
|-|-|-|
|$A,a$|Action|动作,$A$ 表示动作集合,$a$ 表示动作实例|
|DP|Dynamic Programming|动态规划|
|$G$|Gain, Returns|回报,是收益的累积|
|MC|Mento Carlo|蒙特卡洛方法|
|MP|Markov Property|马尔可夫性质|
|MRP|Markov Reward Process|马尔科夫奖励过程|
|MDP|Markov Decision Process|马尔科夫决策过程|
||On-Policy|同轨策略,同策略,在线策略|
||Off-Policy|离轨策略,异策略,离线策略|
|$\pi$|Policy|策略|
|$\pi_*$|Best Policy|最优策略|
|$Q,q$|Action Value|动作价值函数|
|$Q_\pi,q_\pi$|Action Value under Policy|策略 $\pi$ 下的动作价值函数|
|$Q_*,q_*$|Optimal Action Value|最优(策略)动作价值函数|
|$R,r$|Reward|收益或奖励|
|$S,s$|State|状态,$S$ 表示状态集合,$s$ 表示状态实例|
|$V,v$|State Value|状态价值函数|
|$V_\pi,v_\pi$|State Value under Policy|策略 $\pi$ 下的状态价值函数|
|$V_*,v_*$|Optimal State Value|最优(策略)状态价值函数|
大小表示集合,小写表示实例。

Просмотреть файл

@ -0,0 +1,2 @@
# 第 10 章 冰面行走问题

Просмотреть файл

@ -0,0 +1,230 @@
## 10.1 冰面行走问题
在 10.2 节中,我们重温了第六章和第七章中的相关内容。在第八章中,引入了策略问题,并使用贝尔曼期望方程来评估策略的价值。在本章中后续的内容,就是要学习在我模型的情况下,如何使用蒙特卡洛法来评估策略的价值,即估算状态价值函数 $V_\pi$ 和动作价值函数 $Q_\pi$。
### 10.1.1 提出问题
冰面行走问题是 OpenAI Gym 中的一个强化学习环境。在 Gym 文档 https://www.gymlibrary.ml/ 中,提供了以下一些环境:
- Atari - 雅达利游戏
包含一些简单的射击、追逐、对战游戏,玩家通过学习环境变化来控制物体执行各种动作以赢得游戏。
- MuJoCo- 多关节运动控制。
Multi-Joint dynamics with Contact具有接触的多关节动力学。它是一个物理模拟器用于促进机器人、生物力学、图形和动画以及其他需要快速准确模拟的领域的研究和开发。
- Toy Text - 玩具文本环境
都是使用本机Python库如StringIO创建的。这些环境设计得非常简单具有很少的离散状态和动作空间因此易于学习适用于强化学习算法的调试实现。所有环境都可以通过每个环境文档中指定的参数进行配置。
- Classic Control - 经典控制环境
提供了五种典型的控制环境,所有这些环境在给定范围内的初始状态都是随机的。这组环境更容易通过策略来解决,可以通过每个环境文档中指定的参数进行配置。
- Box2D - 二维仿真控制游戏
这些环境都涉及基于物理控制的玩具游戏,使用基于 box2d 的物理和基于 PyGame 的渲染。这些环境是奥列格·克利莫夫在早期创造的,从那时起就成为了流行的玩具基准。所有环境都可以通过每个环境文档中指定的参数进行配置。
我们要解决的问题叫做 Frozen Lake直译为冰冻湖面即冰面行走归属于 Toy Text 环境,如图 10.1.1 所示:
<center>
<img src='./img/FrozenLake-Start.png'>
图 10.1.1 FrozenLake问题示意图
(左图:图形化界面;右图:字符化界面)
</center>
为什么叫做 Toy Text 呢?看图 10.1.1 的右图就可以明白了,这些问题是可以用字符化界面来展示说明并模拟推演的,其中第一行第二个字符 F 底色是灰色的,表示智能体当前的位置。字符代表的含义解释如下。
在一个 4x4 网格的冰冻湖面上,存在 4 种状态区域:
- S - Start 起点,一个,左上角;
- F - Frozen 冰面,十个,图中浅蓝色网格;
- H - Hole 冰洞,四个,图中深蓝色网格;
- G - Goal 目的地,一个,右下角。
环境说明:
- 目标:让图中的小人(智能体)学会从起点尽量安全地走到终点拿到礼物,不要掉入冰洞。
- 分幕:一旦掉入冰洞,本幕结束,回到起点重新开始。拿到礼物也结束一幕。
- 奖励:只有到达目的地才有 +1 的奖励,其它任何状态、动作都是 0 分的奖励(相当于没有奖励)。
- 动作:上下左右四个动作。出界时回到动作出发点。
另外,在冰面行走时有两种难度选项:
- 简单模式:智能体会沿着指定的行走方向顺利地到达下一格。如图 10.1.2 左图所示。
- 困难模式:由于冰面很滑,智能体有 1/3 的概率偏向出发点与动作方向左侧的格子1/3 的概率偏向右侧的格子,只有 1/3 的概率会沿着指定的行走方向到达下一格。如图 10.1.2 右图所示,实心箭头表示行走方向,虚线箭头表示有可能到达的方向。
<center>
<img src="./img/FrozenLake-2Cases.png">
图 10.1.2 冰面行走问题的两种难度模式
(左图:简单模式;右图:困难模式)
</center>
比如图 10.1.2 中智能体所在的位置为 $s_6$,如果向上方走,在简单模式下,会直接走到 $s_2$;在困难模式下有可能会滑倒左侧 $s_5$ 或右侧 $s_7$ 冰洞中,只有 1/3 的概率可以到达 $s_2$。
### 10.1.2 有模型与无模型
有的读者可能会觉得这个问题很简单,只需要用代码判断下一个格子是不是冰洞,如果是的话就绕开,然后让智能体尽量向右下方移动,就可以完成任务了。为什么要用强化学习?
如同所有的机器学习一样,强化学习就是要避免使用“程序”来控制智能体的行为,而是让它真正具有一些智能。
曾经有人讨论电饭锅是不是智能体,因为它可以自动地把米饭做熟。电饭锅里面只是一个 Program程序不是 Intelligence智能比如它不会知道锅里放了多少米和多少水。它只有一个温度传感器和一个定时器当在一个温度下持续了 30 分钟后,它就会认为米饭熟了。
而强化学习的解题思路是计算状态价值函数和动作价值函数,让智能体“学习”到在什么状态下要执行什么动作,才能获利最大。冰面行走问题是一个标准的马尔可夫决策过程 MDP 的场景如第八章所学MDP 下有两种价值函数:状态价值函数和动作价值函数。
$$
v_\pi(s) = \mathbb E[G_t \mid S_t=s]=\mathbb E \Big[\sum_{k=0}^T \gamma^k R_{t+k+1} \mid S_t=s\Big] \tag{10.1.1}
$$
$$
q_\pi(s,a) = \mathbb E[G_t \mid S_t=s, A_t=a]=\mathbb E \Big[\sum_{k=0}^T \gamma^k R_{t+k+1} \mid S_t=s,A_t=a\Big] \tag{10.1.2}
$$
10.1.1和式10.1.2)与 8.4 节中的关于 $v_\pi,q_\pi$ 的表达式有所不同,这里没有写出后续的含有模型信息的公式,即状态转移概率,这是因为在实际问题中,有很多场景是得不到状态转移概率的,无法做精确计算,只能使用蒙特卡洛方法来做“近似”。
有模型的问题,即 model-based如果模型已知则可以用动态规划来解决叫做 Planning。因为是迭代计算所以不算是数学解析解只能叫做高精度的迭代解。
无模型的问题,即 model-free用本章学习的蒙特卡洛法来解决它不试图得到环境模型的数学描述而是通过大量样本的奖励值得到“平均”值叫做近似解。
高方差、低偏差。。。。。。
### 10.1.3 用动态规划法得到近似解
#### 冰面行走问题的两种模式
这就又回到以前的问题:如果只能近似的话,那么以什么为基准来做近似呢?如何判断近似得到位了呢?
对于 Gym 中的这个冰面行走问题,其实它的背后是有一个环境模型的,而 Gym 中的其它很多环境,比如雅达利游戏,就属于没有环境模型的,或者是环境模型过于复杂无法描述的,比如 21 点纸牌游戏。所以,我们可以先用第八章学习的贝尔曼期望方程和第九章学习的贝尔曼最优方程,来得到冰面行走问题的“迭代”解,然后再用蒙特卡洛方法来“近似”,得到一些基本的概念,然后再去解决无模型的问题。
为了使用动态规划法,需要先查看环境提供的数据结构。
【代码位置】MC_103_FrozenLake_DP.py
```python
env = gym.make('FrozenLake-v1', map_name = "4x4", is_slippery=False) # 设置环境为“简单”
print(env.P) # env.unwrapped.P
```
得到如下结果:
```
{
0: {0: [(1.0, 0, 0.0, False)], 1: [(1.0, 4, 0.0, False)], 2: [(1.0, 1, 0.0, False)], 3: [(1.0, 0, 0.0, False)]},
1: {0: [(1.0, 0, 0.0, False)], 1: [(1.0, 5, 0.0, True)], 2: [(1.0, 2, 0.0, False)], 3: [(1.0, 1, 0.0, False)]},
......
14: {0: [(1.0, 13, 0.0, False)], 1: [(1.0, 14, 0.0, False)], 2: [(1.0, 15, 1.0, True)], 3: [(1.0, 10, 0.0, False)]},
15: {0: [(1.0, 15, 0, True)], 1: [(1.0, 15, 0, True)], 2: [(1.0, 15, 0, True)], 3: [(1.0, 15, 0, True)]}
}
```
- 首先这是一个字典形式的数据结构key 是状态序号0-15
- 其次,每个状态序号下面的 value 是另外一个字典,其 key 是动作序号0-3
- 第三级,动作序号下面的 value 是一个列表,里面只有一项用 tuple 表示的数据比如1.0, 0, 0.0, False其含义是
- 1.0:动作发生后的转移概率;
- 0转移到的目标状态序号
- 0.0:得到的奖励值;
- False目标状态是否为终止状态。
如果设置环境为“困难”,即把参数 is_slippery 设置为 True
```python
env = gym.make('FrozenLake-v1', map_name = "4x4", is_slippery=True) # 设置环境为“困难”
print(env.P) # env.unwrapped.P
```
得到如下结果:
```
{
0: { # 状态 0
0: [(0.333, 0, 0.0, False), (0.333, 0, 0.0, False), (0.333, 4, 0.0, False)], # 动作0
1: [(0.333, 0, 0.0, False), (0.333, 4, 0.0, False), (0.333, 1, 0.0, False)], # 动作1
2: [(0.333, 4, 0.0, False), (0.333, 1, 0.0, False), (0.333, 0, 0.0, False)], # 动作2
3: [(0.333, 1, 0.0, False), (0.333, 0, 0.0, False), (0.333, 0, 0.0, False)] # 动作3
},
......
14: { # 状态 14
0: [(0.333, 10, 0.0, False), (0.333, 13, 0.0, False), (0.333, 14, 0.0, False)], # 动作 0
1: [(0.333, 13, 0.0, False), (0.333, 14, 0.0, False), (0.333, 15, 1.0, True)], # 动作 1
2: [(0.333, 14, 0.0, False), (0.333, 15, 1.0, True), (0.333, 10, 0.0, False)], # 动作 2
3: [(0.333, 15, 1.0, True), (0.333, 10, 0.0, False), (0.333, 13, 0.0, False)] # 动作 3
},
15: { # 状态 15
0: [(1.0, 15, 0, True)], # 动作0 (转移概率, 下个状态, 奖励, 是否终止)
1: [(1.0, 15, 0, True)], # 动作1
2: [(1.0, 15, 0, True)], # 动作2
3: [(1.0, 15, 0, True)] # 动作3
}
}
```
仍然是一个二级字典结构,但是在第三级数据中加入了状态转移的描述,比如状态 0 下的动作 2 下列表由三组 tuple 组成,表示在状态 0 执行动作 2向右移动
- (0.333, 4, 0.0, False) ,移动方向的右侧,**滑到**状态 4 的概率为 0.333,奖励 0.0,非终止状态;
- (0.333, 1, 0.0, False) ,移动方向的前方,**达到**状态 1 的概率为 0.333,奖励 0.0,非终止状态;
- (0.333, 0, 0.0, False) ,移动方向的左侧,**回到**状态 0 的概率为 0.333,奖励 0.0,非终止状态。
注意我们使用了三种不同的词汇来描述状态转移:
- **滑到**:表示与动作预期的方向不同;
- **达到**:表示与动作预期的方向相同;
- **回到**:本来是应该向上方滑动,但是出界了,所以回到了当前的位置。
状态 15 的数据比较简单,和环境为“简单”时相同,因为该状态是终止状态,无论任何动作都会回到原地。但是细心的读者会发现,该状态中并没有得到奖励为 +1 的信息,这是为什么?
因为我们前面讲“奖励”的时候说过,正确的奖励定义是面向转移过程的。在这个问题中,只有从状态 14 才有可能到达状态 15因此我们可以再看一下状态 14 的数据:
- 动作 0因为是向左移动所以没有任何可能达到状态 15
- 动作 1向下移动有可能向左**滑到**状态 15奖励 1.0,列表第三项;
- 动作 2向右移动有可能直接**达到**状态 15奖励 1.0,列表第二项;
- 动作 3向上移动有可能向右**滑到**状态 15奖励 1.0,列表第一项。
#### 冰面行走问题简单模式下的随机策略的迭代解
有了上述对数据的理解,把第九章的代码稍微改动一下,就可以完成对 $v_\pi,q_\pi$ 的迭代计算了。
【代码位置】Algorithm.Algo_PolicyValueFunction.py, MC_103_FrozenLake_DP.py
$V_\pi$ 结果如下:
```
====================is slippery=False====================
迭代次数 = 21
---------------随机策略下的状态价值函数 V---------------
[[0.014 0.011 0.021 0.01 ]
[0.016 0. 0.041 0. ]
[0.035 0.088 0.142 0. ]
[0. 0.176 0.439 0. ]]
```
$Q_\pi$ 结果如下:
```
---------------随机策略下的动作价值函数 Q---------------
[[0.014 0.016 0.011 0.014] +-----+-----+-----+-----+
[0.014 0. 0.021 0.011] | | | | |
[0.011 0.041 0.01 0.021] | ┼ | ┼─►| ┼ |◄─┼ |
[0.021 0. 0.01 0.01 ] | ▼ | | ▼ | |
[0.016 0.035 0. 0.014] +-----+-----+-----+-----+
[0. 0. 0. 0. ] | | | | |
[0. 0.142 0. 0.021] | ┼ | ┼ | ┼ | ┼ |
[0. 0. 0. 0. ] | ▼ | | ▼ | |
[0.035 0. 0.088 0.016] +-----+-----+-----+-----+
[0.035 0.176 0.142 0. ] | | | | |
[0.088 0.439 0. 0.041] | ┼─►| ┼ | ┼ | ┼ |
[0. 0. 0. 0. ] | | ▼ | ▼ | |
[0. 0. 0. 0. ] +-----+-----+-----+-----+
[0. 0.176 0.439 0.088] | | | | |
[0.176 0.439 1. 0.142] | ┼ | ┼─►| ┼─►| ┼ |
[0. 0. 0. 0. ]] +-----+-----+-----+-----+
```
从 $Q$ 函数结果以及抽取出的策略(方格内的箭头)来看,这个结果是非常正确的,可以引导智能体从起点安全地到达终点。
有了上述基准后,我们就可以去研究最佳的算法了。

Просмотреть файл

@ -0,0 +1,248 @@
## 10.2 蒙特卡洛方法
### 10.2.1 提出问题
前面的章节以贝尔曼方程 $\to$ 贝尔曼期望方程 $\to$ 贝尔曼最优方程为线索,用动态规划思想解决强化学习问题,其基础是我们知道了环境的状态转移概率 $P$ 和奖励模型 $R$未来变成了可预测的。这种情况称作是有模型model-based问题。
但是在现实生活中,很多情况下无法给出一个准确的 $P$,也不知道奖励函数 $R$,原因可能是状态太多,或者转移太复杂,以至于无法建模或者计算量巨大。
所以问题是:**当无模型时model-free无法使用动态规划如何计算状态价值和动作价值以及策略**
这就要依赖在本章中将要学习的**蒙特卡洛方法**Monte Carlo Method它是强化学习的一种重要手段。其实在“蒙提霍尔问题”和“安全驾驶问题”中已经使用了蒙特卡洛的思想来解决问题了。在后续的章节中要充分地讨论蒙特卡洛预测和蒙特卡洛控制都是以本节的内容为理论基础的。
接下来我们系统地了解一下这种方法。
### 10.2.2 应用场景一
#### 古人如何计算圆周率?
假设我们已经知道了圆面积的方程为 $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,艰难地把圆周率推到了小数点后四位。
- 公元 480 年前后
祖冲之继续使用割圆术计算12,288形的边长将圆周率推到了小数点后六位可惜的是由于文献的失传祖冲之的计算方法我们现在已经不得而知了。
- 十五世纪
随着近代数学的发展,数学家韦达、罗门、科伊伦、司乃耳、格林伯格通过割圆术陆续将圆周率推到了小数点后 39 位,这个精度是什么概念呢,如果我们通过小数点后 39 位的圆周率计算一个一个可观察宇宙大小的圆,计算的误差仅仅只有一个氢原子大小。
- 十六世纪到十七世纪
人们发现了一种新的圆周率计算方法——无穷级数法让计算圆周率的工作变得更加快速。无穷级数是一组无穷数列的和数学家梅钦通过无穷级数将圆周率推算到小数点后100位在很短的时间里人们通过梅钦类公式反复打破了新的圆周率记录。
- 十八世纪
法国数学家布丰提出了随机投针法,即利用概率统计的方法来计算圆周率的值,也就是著名的投针实验。布丰在地板上画出若干平行的直线,再将一根根短于平行直线距离的针撒到地板上,通过统计针的总数和与直线相交的针的个数,从而计算圆周率。
这种算法虽然虽然没有打破圆周率的记录,但这种将几何与概率结合起来的思想催生了蒙特卡洛算法,也让人工智能成为了可能。
#### 用蒙特卡洛方法估算圆周率
先给出一个试验方法,在圆形外面套一个正方形,如图 10.2.1 所示。
<center>
<img src="./img/CircleSquare.png">
图 10.2.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$。
如果是古人的话,可能会在图 10.2.1 上平铺一层细沙,然后称出在圆内的细沙的重量,再称出方形内细沙的重量,以此来模拟面积计算。现在我们借助计算机的力量,在图中随机打点,然后统计出圆内的点的数量和方形内点的数量,也同样可以模拟计算出面积。
#### 代码实现
【代码位置】MC_101_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()
```
为了避免随机性,可以使用相同的点数做多次试验并求平均,或者使用不同的点数做多次试验求平均。在此我们使用了第二种方法,绘制出图 10.2.2 所示的结果。
<center>
<img src="./img/CirclePi-4.png">
图 10.2.2 随机投点估算圆周率
</center>
在本次试验中,随机点数量与估算出的 $\pi$ 值对应关系如下:
- 100个点3.08
- 200个点3.02
- 500个点3.176
- 1000个点3.156
- 平均3.108
人们通常以为点数越多,估算出的值越准确。从趋势上来看确实如此,但也不能避免单次试验的结果都带来的偏差。因此,在图 10.2.3 的试验中,我们使用 [1000, 1100, 1200, ..., 20000] 个点来做估算,并且每个点数上都重复了 100 次取平均,最后得到了误差的变化图。
【代码位置】MC_101_CirclePi_2.py
<center>
<img src="./img/CirclePi-error.png">
图 10.2.3 随机点计算圆周率的误差变化
</center>
从横坐标上看投点数量越多偏差和方差都会随之减小。从最终的平均结果来看3.1417 已经可以代表这种方法的精度了。当然,读者还可以用更多的点数来验证,不过单次试验的结果总会令人失望的。
对于估算不规则图形面积的问题,用这种方法同样可以解决。
### 10.2.3 应用场景二
#### 复杂函数的定积分如何计算?
对于复杂函数,如何计算其定积分?
在图 10.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$
<center>
<img src="./img/Integral.png">
图 10.2.4 计算函数积分
</center>
前两个简单函数,可以得知其定积分结果为:
$$
\int_{0.2}^1 \frac{1}{x^2} dx = 4, \quad \int_0^{\pi} \sin(x) dx = 2
$$
那么对第三个函数做积分运算就比较麻烦了。
#### 用蒙特卡洛方法估算定积分值
从问题一中,我们学到了随机投点的方法,在计算定积分时仍然可以使用,但是由于函数值域区间的变化,写出代码来有些复杂(但还是可以完成的)。所以这一次介绍一个不同的方法,见图 10.2.5。
<center>
<img src="./img/Integral-2.png">
图 10.2.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{10.2.1}
$$
如果 $N$ 足够大的话,$S$ 的值会非常接近于数学积分值。
其理论基础是,随机变量 $X$ 的数学期望为:
$$
\mathbb E[X]=\int_{-\infty}^{\infty} pdf(x)·x \ dx \tag{10.2.2}
$$
其中 $pdf(x)$ 是概率分布函数,变成离散型公式就是:
$$
S = \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{pdf(x_i)} \tag{10.2.3}
$$
如果我们使用均匀分布,则 $pdf(X_i)=\frac{1}{b-a}$是一个常量所以式10.2.3可以变形为式10.2.1)。
#### 代码实现
【代码位置】MC_101_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
# 求积分 式10.2.1
def integral(f, a, b, n):
v = 0
repeat = 10
for i in range(repeat):
x = np.random.uniform(a, b, size=(n, 1)) # 在[a,b]区间上均匀随机取点
y = f(x)
v += np.sum(y) / n * (b-a) # 式(10.2.2)
return v/repeat
```
n=100 时,最后得到的结果是:
```
S1 = 3.9978943535545155
S2 = 1.9983227898197693
S3 = -0.1512574314084149
```
近似为 4.00, 2.00, -0.15。
### 10.2.4 蒙特卡洛方法
蒙特卡洛Monte Carlo, MC方法也称统计模拟方法是一类随机算法的统称是在 1940 年后由于科学技术的发展和电子计算机的发明,而提出的一种以概率统计理论为指导的数值计算方法。是指使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。
20 世纪 40 年代,在冯·诺伊曼、斯塔尼斯拉夫·乌拉姆和尼古拉斯·梅特罗波利斯在洛斯阿拉莫斯国家实验室为核武器计划工作时,发明了蒙特卡罗方法。因为乌拉姆的叔叔经常在摩纳哥的蒙特卡洛赌场输钱得名,而蒙特卡罗方法正是以概率为基础的方法。与它对应的是确定性算法。
蒙特卡洛方法是一种近似推断的方法,通过采集大量样本的方法来求解期望、均值、面积、积分等问题。蒙特卡洛对某一种分布的采样方法有直接采样、接受拒绝采样与重要性采样三种,直接采样最简单,但是需要已知累积分布的形式;接受拒绝采样与重要性采样适用于原分布未知的情况,这两种方法都是给出一个提议分布,不同的是接受拒绝采样对不满足原分布的样本予以拒绝,而重要性采样则是给予每个样本不同的权重,大家可以根据不同的场景使用这三种方法中的一种进行采样。
读者不要单纯地认为点越多越精确,这是有一个前提的。总体而言,不精确是因为我们的点还没有到某个量级。一般情况下,蒙特卡罗算法的特点是,采样越多,越近似最优解,而永远不是最优解。
其优点比较明显,尤其对于具有统计性质的问题可以直接进行解决,对于连续性的问题也不必进行离散化处理。其缺点也很明显,对于确定性问题转化成随机性问题做的估值处理,丧失精确性,得到一个接近准确的值也不太容易。
蒙特卡洛方法可以解决的问题可以粗略地分成两类:
- 所求解的问题本身具有内在的随机性,借助计算机的运算能力可以直接**模拟随机的过程**
例如在核物理研究中,分析中子在反应堆中的传输过程。中子与原子核作用受到量子力学规律的制约,人们只能知道它们相互作用发生的概率,却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向。
科学家依据其概率进行随机抽样得到裂变位置、速度和方向,这样模拟大量中子的行为后,经过统计就能获得中子传输的范围,作为反应堆设计的依据。
- 所求解问题可以转化为某种**随机分布的特征数**
比如随机事件出现的概率,或者随机变量的期望值。通过随机抽样的方法,以随机事件出现的频率估计其概率,或者以抽样的数字特征估算随机变量的数字特征,并将其作为问题的解这种方法多用于求解复杂的多维积分问题。
我们在上面举的例子是一个一维的积分问题,积分结果是曲线下的面积;当然可以再增加一维随机采样,变成二维积分问题,结果为曲面下的体积,比一维积分更有实际作用。
蒙特卡洛在强化学习中应用的核心主要包含以下几点:
- 是直接从经验轨迹当中直接进行学习。
- 是一种无模型model-free方法即没有 MDP 的转移概率以及奖励的先验知识。
- 从完整的经验轨迹中学习不使用自举bootstrapping方法这与前几章学习的动态规划方法不同。
- 简单地使用这个思想:期望值 = 平均回报。即采样后得到回报并累积,最后除以采样次数等到平均值。
- 缺点在于只能应用于一定有终结点的分幕的 MDP 过程。

Просмотреть файл

@ -0,0 +1,349 @@
## 10.3 价值评估
在 6.3 节中,已经学习过如何使用最朴素的蒙特卡洛方法来评估价值函数。本小节学习使用蒙特卡洛法估算 MRP 下的价值函数 $V$,将要接触到一种算法:首次访问法。需要说明的是,这种算法同样适用于估算 MDP 下的价值函数 $V_\pi, Q_\pi$。
### 10.3.1 提出问题
#### 安全驾驶问题中的状态价值函数估算方法
首先要回忆一下,在 MRP 中,因为没有牵涉到动作问题,所以只有状态价值函数的概念。
在第六章中,针对 MRP 问题,我们曾经使用蒙特卡洛方法估算过“安全驾驶问题”中的状态价值函数。当时给读者留下的印象可能是:
- 这种方法怎么运行起来那么慢?
- 因为要采样很多次才能基本稳定下来,否则两次试验之间的差别特别大,那么到底采样多少次才算合理?
- 这种方法怎么能保证准确?有没有基准可以比较?
- 最终得到的结果的可信度有多高?
在 6.3 节中,我们根据回报 $G$ 的定义:
$$
G_t = R_{t+1}+\gamma R_{t+2}+\gamma^2 R_{t+3}+ \cdots +\gamma^{T-t-1} R_{T}
\tag{10.3.1}
$$
以及价值函数 $V$ 的定义:
$$
v_t(s) = \mathbb E [G_t | S_t = s]
\tag{10.3.2}
$$
然后使用最朴素的蒙特卡洛采样法初步计算出了安全驾驶问题的各个状态的价值函数。简述其过程如下:
1. 使用蒙特卡洛采样,每次采样都需要指定一个初始状态 $s_0$,然后开始在幕内循环;
2. 得到一个 $R$ 后,累积计算 $G = G + \gamma^{t-1}R$,直到幕内循环结束(到 $T$得到式10.3.1)的计算结果;
3. 进行下一次采样,$G$ 值清零后再用 step 2 的循环重新计算 $G$ 值并累积;
4. 最后求 $G$ 的平均(数学期望)得到 $v_t(s_0)$即式10.3.2)。
这种方法的特点是:
- 易理解,完全按照定义来计算;
- 省空间,不需要记录中间状态,来了一个 $R$ 后立刻计算,然后扔掉,只保留 $G$ 值;
- 随着幕数的增加,会更逼近真实值;
- 每次多幕循环只计算一个起始状态的价值;
- 遍历状态并在环境控制下重复多次,速度慢;
- 浪费了中间的采样结果。
#### 充分利用样本数据
聪明的读者可能会发现一个问题:在得到一幕的采样后,如果不“从头”开始,而是从第二个、第三个状态开始计算,是不是就能在一次采样中就可以得到很多状态的 $G$ 值呢?
<center>
<img src='./img/MC-1.png'>
图 10.3.1 重复利用样本来计算不同时刻/状态的回报值
</center>
如图 10.3.1 所示,从 $S_0$ 开始一幕的采样,到 $T$ 为止结束,得到 $R_1,\cdots,R_T$ 的序列后:
- 第一行:固然可以从 $R_1$ 开始计算出 $G_0$
- 第二行:但是如果从 $R_2$ 开始,不就能计算出 $G_1$ 了吗?
- 第三行:同理,还可以计算出这一采样序列中的任意的 $G_t$ 出来。
这样的话利用一次采样结果可以计算出很多状态的 $G$ 值,会大幅提高算法的效率。当然,这就需要把每一幕的采样结果记录下来,计算完所有可能的 $G$ 后,扔掉这一幕记录,再进行下一幕的采样。
#### 用贝尔曼方程得到解析解
在第六章时,我们还没有学习贝尔曼方程(第七章的内容),所以当时使用了蒙特卡洛法依靠大量的采样来计算 $G$ 值。由于安全驾驶问题实际上是有一个模型的所以可以利用贝尔曼方程得到该问题的“近似真实解”并以此为基准Ground Truth来衡量算法的性能。而性能又包括两个方面1. 速度2. 准确度。需要在这二者之间寻找平衡。
下面是用贝尔曼方程的矩阵法计算出来的状态价值函数值 $V$(精确到小数点后两位)。
【代码位置】MC_102_SafetyDrive_DataModel.py
```
状态价值函数计算结果(数组) : [ 1.03 1.72 2.72 3.02 -5.17 -6.73 6. -2.37 -1. 5. 0. ]
Start: 1.03
Normal: 1.72
Pedestrians: 2.72
DownSpeed: 3.02
ExceedSpeed:-5.17
RedLight: -6.73
LowSpeed: 6.0
MobilePhone:-2.37
Crash: -1.0
Goal: 5.0
End: 0.0
```
但是我们也讲到过,当状态数量很多时,用矩阵法计算的复杂度很高,而且很多问题没有模型,因此也不能使用贝尔曼方程(动态规划)。
#### 衡量算法准确度 RMSE
RMSE - Root Mean Square Error均方根误差用于比较观测值与真实值之间的偏差通常不是比较两个标量而是比较两个向量数组。比如可以用贝尔曼方程的矩阵法计算出来的状态价值数组做为基准值 $Y$,用蒙特卡洛法估算出来的数组值 $X$与基准值 $Y$ 之间计算 RMSE。
$$
RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^n (x_i-y_i)^2} \tag{10.3.3}
$$
其中:
- $n$ 为数组元素个数;
- $x_i$ 为观测值;
- $y_i$ 为真实值。
代码实现为:
```python
def RMSE(x, y):
err = np.sqrt(np.sum(np.square(x - y))/y.shape[0])
return err
```
比如,假设真实值为 $y=[1,2,3]$,观测值为 $x=[1.1,2,2.9]$,则:$RMSE=\sqrt{[(1.1-1)^2+(2-2)^2+(2.9-3)^2]/3} \approx 0.08$。
在做向量比较时,有一些资料上把 $n$ 解释为观测次数其实是不对的。另外在式10.3.3)中的 $x,y$ 是可以对调的,不影响计算结果。
注意,不同的问题之间的 RMSE 是不可比的,比如在问题 A 中,理想的 RMSE 应该小于 0.01,而在问题 B 中,只要小于 0.1 就达到要求了。最简单的做法是先用全 0 值与基准值作比较,计算出 RMSE_0然后把该值除以 50 得到 RMSE_0_50只要后续的算法达到 RMSE_0_50 的级别,就算符合要求了。
用 $y=[1,2,3]$ 举例来说RMSE_0=$\sqrt{[(1-0)^2+(2-0)^2+(3-0)^2]/3}/50 \approx 0.0432$,只要后续算法的误差为这个级别就行了。
如果用安全驾驶问题的基准值 y=[1.03, 1.72, 2.72, 3.02, -5.17, -6.73, 6, -2.37, -1, 5, 0] 来做例子的话RMSE_0_50 = 0.076。
### 10.3.2 首次访问法
#### 什么是首次访问法
在图 10.3.1 中,$s_0,s_1,\cdots,s_t$ 代表的是不同时刻的状态,而不是不同的状态。说起来比较绕,我们举个例子。在安全驾驶问题中,一个序列可能是这样的:出发$s_0 \to$正常行驶$s_1 \to$闹市减速$s_2 \to$礼让行人$s_3\to$正常行驶$s_4\to \cdots$
读者可以看到,正常行驶状态出现了两次,一次在 $s_1$,一次在 $s_4$,对这个状态来说,$s_1$ 就叫做首次访问。在这两个状态可以计算出 $G_1$ 和 $G_4$ 来,那么最后 $V_{正常行驶}$ 的值是 $G_1$ 呢还是 $G_4$ 呢?
更一般地在一个采样序列中有如下状态序列数据Episode=$[s_a, s_b, s_c, s_a, s_d, s_b]$, $\tau=6$。
- 首次访问法在倒序遍历到最后那个 $s_b$ 时,$t=\tau-1=5$,先检查 $s_b$ 是否在 Episode[0:5] 中;在本例中 Episode[0:5] 含有 $s_b$,所以就不会把这个 $t$ 上的 $G_t$ 值记录到 $G$ 的累积值中。
- $t=4$$s_d$ 在 Episode[0:4] 中不存在,叫做首次访问,记录其 $G_t$ 值;
- $t=3$$s_a$ 在 Episode[0:3] 中存在,不记录其 $G_t$ 值;
......
在这个例子中,最终会计算的是 $s_a, s_b, s_c, s_d$ 四个状态的 $G$ 值。
马尔科夫链的环状结构会使得这种情况非常常见,显然第一个 $s_a$ 和 第四个 $s_a$ 在时域 $t$ 上应该具有不同的价值。但是从全局/静态角度看它们又应该具有相同的价值。Sutton 在研究了这个问题后,认为根据大数定律,首次访问法的方差以 $1/\sqrt{n}$ 的速度收敛。
#### 算法描述
【算法 10.3】首次访问型蒙特卡洛法。
下面的伪代码中,$\leftarrow$ 表示赋值,$\Leftarrow$ 表示追加列表。
---
输入:起始状态 $s,\gamma$, Episodes
初始化:$G(S) \leftarrow 0, N(S) \leftarrow 0$
多幕 Episodes 循环:
  列表置空 $Episode = [\ ] $ 用于存储序列数据 $(s,r)$
  幕内循环直到终止状态:
    从 $s$ 根据环境模型得到 $s',r$ 以及是否终止的标志
    $Episode \Leftarrow (s',r)$
    $s \leftarrow s'$
  $G_t \leftarrow 0$
  对 $Episode$ 从后向前倒序遍历, $t=\tau-1,\tau-2,...,0$
    取出 $s_t,r_t$
    $G_t \leftarrow \gamma G_t+r_t$
    如果 $s_t$ 不在 $s_0,s_1,\cdots,s_{t-1}$ 中,即首次访问:
      $G(s_t) \leftarrow G(s_t)+G_t$
      $N(s_t) \leftarrow N(s_t)+1$
$V(S) \leftarrow G(S) / N(S)$
输出:$V(S)$
---
#### 倒序遍历方法
图 10.3.2 是倒序遍历方法的说明。
<center>
<img src="./img/MC-2.png">
图 10.3.2 首次访问法中计算 G 值的方法
</center>
把这张图与图 6.3.1 做比较,就可以发现橙色部分的“计算顺序”是不同的。其中,状态下标 $_S,_N,_L,_G,_E$ 等是安全驾驶问题中对应的状态的缩写,在这里无关紧要,但它们是一个真实的序列。
“倒序遍历”与“首次访问”,这两个顺序是正好相反的,请读者仔细甄别其具体含义。
#### 算法实现
【代码位置】Algorithm.Algo_MC_Value_Evaluation.py
```Python
# MC1-FirstVisit 首次访问法
def MC_FirstVisit(dataModel, start_state, episodes, gamma):
Value = np.zeros(dataModel.nS) # G 的总和
Count = np.zeros(dataModel.nS) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
TrajectoryState = [] # 一幕内的状态序列
TrajectoryReward = [] # 一幕内的奖励序列
s = start_state
is_end = False
while (is_end is False): # 幕内循环
next_s, r, is_end = dataModel.step(s) # 从环境获得下一个状态和奖励
TrajectoryState.append(s.value)
TrajectoryReward.append(r)
s = next_s
assert(len(TrajectoryState) == len(TrajectoryReward))
num_step = len(TrajectoryState)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s = TrajectoryState[t]
r = TrajectoryReward[t]
G = gamma * G + r
if not (s in TrajectoryState[0:t]):# 是否首次访问
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
return Value / Count # 求均值
```
与每次访问法代码不同的是,上述代码使用了两个 Trajectory 来分别记录状态序列和奖励序列(当然这二者是同步的),便于在判断是否首次访问时简单易行。
#### 运行结果
【代码位置】MC_102_FirstVisit.py
首次访问法的运行结果:
```
gamma = 1.0
状态价值函数计算结果(数组) : [ 1.03 1.71 2.73 2.83 -5.2 -6.71 6. -2.28 -1. 5. 0. ]
Start: 1.03
Normal: 1.71
Pedestrians: 2.73
DownSpeed: 2.83
ExceedSpeed:-5.2
RedLight: -6.71
LowSpeed: 6.0
MobilePhone:-2.28
Crash: -1.0
Goal: 5.0
End: 0.0
耗时 : 4.39701247215271
误差 = 0.06336656060380948
```
#### 几种算法的运行结果比较
表 10.3.1 $\gamma=1$ 时的状态值计算结果比较
|状态|矩阵法<br>(第七章)|顺序计算法<br>(第六章)|首次访问法<br>(本章)|
|-|:-:|:-:|:-:|:-:|
|出发 Start| 1.03|1.11|1.03|
|正常行驶 Normal| 1.72|1.52|1.71|
|礼让行人 Pedestrians| 2.72|2.70|2.73|
|闹市减速 DownSpeed| 3.02|2.79|2.83|
|超速行驶 ExceedSpeed| -5.17|-5.19|-5.2|
|路口闯灯 RedLight| -6.73|-6.71|-6.71|
|小区减速 LowSpeed| 6.00|6.00|6.00|
|拨打电话 MobilePhone| -2.37|-2.51|-2.28|
|发生事故 Crash| -1.00|-1.00|-1.00|
|安全抵达 Goal| 5.00|5.00|5.00|
|终止 End| 0.00| 0.00|0.00|
|**误差 RMSE**|-|**0.102**|**0.063**|
|**总采样次数**|-|**50000**|**100000**|
|**耗时**|-|**7.37**|**4.70**|
可以看到改进的算法在速度上和精度上都比原始算法要好。RMSE 的误差值,原始顺序计算法误差为 0.102(大于前面说过的 RMSE_0_50=0.076),改进算法为 0.063(小于 0.076),越小越好。
从性能上看,第六章介绍的顺序计算法对每个状态做了 5000 次采样,一共 10 个有效状态,所以就是 50000 次采样。首次访问法对所有状态(混合)一共做了 20000 次采样,以每幕平均 5 个状态来估算,就有至少 100000 次采样,结果自然会比顺序计算法要准确。
几种算法的耗时值其实并不准确,因为每次运行都有不同的值,但是可以确定的是,首次访问法比顺序计算法性能要好,因为重复利用了采样结果。
### 10.3.3 何时可以停止迭代
在问题规模比较小时,我们当然可以循环尽可能多的幕数,得到更多的采样来计算均值,以逼近期望值。但是,问题规模较大较复杂时,循环会很耗费时间。所以,如何权衡循环时间和结果精度之间的关系,是一个需要研究的问题。
这可能让读者想起了在动态规划迭代中所采用的方法:当前后两次结果的差值的最大值小于一个给定的阈值时,就可以停止迭代。代码是这样的:
```python
if abs(V-V_old).max() < delta:
break
```
- 动态规划方法有个特点:它的每一次迭代,肯定会改变计算值,而且这种改变是由大变小的。当上述代码成立时,并假设不终止迭代,那么下一次的迭代结果的差值,也一定会小于 delta。
- 而蒙特卡洛法也有个特点:它的每一次循环,也肯定会改变计算值,但是改变的幅度是不可预知的。因为采样是随机的,每次采样可能会改变多个值,也可能只改变一个值。
表 10.3.1
|状态|$s_1$|$s_2$|$s_3$|
|-|-|-|-|
|当前 V 值|0.8|0.9|1.0|
|被访问次数|1000|500|100|
假设新的一幕采样使得 $s_1$ 的被访问次数 +1而且 reward=1那么最新的状态值为$v(s_1)=\frac{0.8 \times 1000+1}{1000+1} \approx 0.8002$,那么 $0.8002-0.8=0.0002 < $1e-3可以看到它的变化非常小但是并不能以这一次的结果来判定应该结束迭代了因为这种采样的变化是累积的只有当多次采样后的变化幅度仍然很小时才可以认定是收敛了。
还有一个问题是阈值的设置,在动态规划问题中,将 delta 设置为 1e-4 甚至 1e-6都是可行的。但是在蒙特卡洛方法中1e-2 或 1e-3 才有可能是一个合理的值。
第三个问题是,由于蒙特卡洛方法的不稳定性,某一次的收敛有可能是假象,我们可以设定为必须连续三次检查都收敛,才算是真的收敛。
检查是否收敛部分的代码:
```python
# MC - FirstVisit with EarlyStop 首次访问法,早停
def MC_FirstVisit_EarlyStop(dataModel, start_state, episodes, gamma, checkpoint=100, delta=1e-3):
......
if (episode+1)%checkpoint == 0:
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
V = Value / Count # 求均值
if abs(V-V_old).max() < delta: # 判断是否收敛
less_than_delta_count += 1 # 计数器 +1
if (less_than_delta_count == 3): # 计数器到位
break # 停止迭代
else: # 没有持续收敛
less_than_delta_count = 0 # 计数器复位
V_old = V.copy()
return Value / Count # 求均值
```
在主程序代码中的几个关键值设置:
```python
episodes = 100000 # 最多循环 100000 次
checkpoint = 200 # 每 200 次检查一次收敛情况
delta = 1e-2 # 阈值为 1e-2
```
运行结果如下:
```
循环幕数 = 22400
状态值 = [1.068 1.742 2.804 3.168 -5.131 -6.672 6. -2.366 -1. 5. 0.]
误差 = 0.057937657059503364
```
一共循环了 22400 次,达到了 3 次收敛的条件,最终的误差为 0.058。
在代码中还记录了每次检查时的 V 值,并计算出与基准值的误差,显示在图 10.3.3 中。
<center>
<img src="./img/EarlyStop.png">
图 10.3.3 早停法
</center>
从图中可以看到,误差的变化幅度还是比较大的,但是到了 23000 次以后,估计会很平缓了,读者可以自行测试。

Просмотреть файл

@ -0,0 +1,268 @@
## 10.4 策略评估(预测)$V_\pi$
本节将会针对无模型问题使用蒙特卡洛法做策略评估policy evaluation也叫做预测Predication
### 10.4.1 估算 $v_\pi$ 和 $q_\pi$
在 10.2 节中,我们使用**首次访问法**,估算了 MRP 问题中的状态价值函数 $V$,也在一开始就提到过,这个方法同样适用于估算 MDP 问题中的价值函数 $v_\pi$ 和 $q_\pi$。
与 MRP 不同的是,在 MDP 中由于引入了**策略**和**动作**的概念,所以,在当前状态 $s$ 下要根据策略policy选择动作action执行动作与环境交互得到下一个状态、奖励值等信息。
表 10.4.1 MRP 与 MDP 的比较
||马尔科夫奖励过程 MRP|马尔科夫决策过程 MDP|
|-|-|-|
|目的|价值评估|策略评估|
|方法|根据 $s$ 直接得到 $s'$|根据 $s$ 查询策略 $\pi \to$ 得到动作 $a$ 并执行 $\to$ 获得奖励 $r$ 和 $s'$。|
|伪代码|s',r = env.step(s)|a = policy(s) # 策略决定动作<br> s',r = env.step(s,a)|
|结果|奖励值累积求平均,得到状态价值函数 $V$|奖励值累积求平均,得到状态价值函数 $V_\pi$ 和动作价值函数 $Q_\pi$|
首先看状态价值函数的定义:
$$
v_\pi(s) \doteq \mathbb E [G_t \mid S_t=s] \tag{10.4.1}
$$
这和式10.2.2)没有区别,所以,算法所需要的参数以及过程都和 10.2 节中一样,只不过是运行在一个 MDP 环境下而已。
另外一个问题是:只有 $V_\pi$ 就足够了吗?答案是否定的。
<center>
<img src="./img/fromV2Q.png">
图 10.4.1 状态价值函数 $V_\pi$ 决定动作价值函数 $Q_\pi$
</center>
如图 10.4.1 中,一个 9 个单元的方格世界中,假设已经知道了 $v(s_1)=1.6,v(s_3)=1.3,v(s_5)=1.2,v(s_7)=1.9$,当智能体处于 $s_4$ 的位置时,它一定应该向具有最大价值的 $s_7$ 移动吗?
不一定看动作价值函数的定义式8.4.2
$$
\begin{aligned}
q_\pi(s,a) & \doteq \mathbb E [G_t \mid S_t = s, A_t=a]
\\
&=P^a_{ss'}[R^a_{ss'}+\gamma V_\pi(s')]
\end{aligned}
\tag{10.4.2}
$$
在有模型的情况下,动作价值由三个因子组成:$R, P, V$。虽然 $v(s_7)$ 最大,但是 $R,P$ 都未知的情况下,不能确定就应该向 $s_7$ 移动。假设 $R^a_{s_4,s_7}=R^a_{s_4,s_5}=1, P^a_{s_4,s_7}=0.5, P^a_{s_4,s_5}=1, \gamma=1$,则:
- $q(s_4,a_{RIGHT})=P^a_{s_4,s_5}[R^a_{s_4,s_5}+\gamma v(s_5)]=1 \cdot [1+1\cdot 1.2]=2.2 $
- $q(s_4,a_{DOWN})=P^a_{s_4,s_7}[R^a_{s_4,s_7}+\gamma v(s_7)]=0.5 \cdot [1+1\cdot 1.9]=1.45 $
由于 $2.2 > 1.45$,那么智能体应该选择向右移动。
在无模型的情况下如果利用式10.4.1)估算出了 $v_\pi$ 后,能直接算出 $q_\pi$ 吗?不能!虽然 $R$ 可以通过与环境交互得到,但是因为不知道 $P$,所以无法计算 $q_\pi$。
在这种情况下我们只能通过式10.4.2)的第一行的定义来估算 $q_\pi$,然后才能按照从 $q_\pi$ 中提出的策略 $\pi(s) \doteq \argmax_a q_\pi(s,a)$ 来行动。所以,在本章后续的内容中,都会把估算 $q_\pi$ 作为目标。
所以,在无模型的情况下,计算动作价值 $q_\pi$ 比起计算状态价值 $v_\pi$ 更有用。在本小节中,我们先学习预测 $V_\pi$,下一小节再学习预测 $Q_\pi$。
### 10.4.2 每次访问法估算 $V_\pi$
#### 算法描述
【算法 10.4】每次访问型蒙特卡洛法。下面的伪代码中,$\leftarrow$ 表示赋值,$\Leftarrow$ 表示追加列表。
---
输入:起始状态 $s$,策略 $\pi$,折扣 $\gamma$, 幕数 Episodes
初始化数组:$G(S) \leftarrow 0, N(S) \leftarrow 0$$S$ 为状态空间
多幕 Episodes 循环:
  列表置空 $Episode = [\ ] $ 用于存储序列数据 $(s,r)$
  幕内循环直到终止状态:
    从 $s$ 根据策略 $\pi$ 得到动作 $a$
    执行 $a$ 从环境得到 $s',r$ 以及是否终止的标志
    $Episode \Leftarrow (s,r)$,相当于是 $(s_t,r_{t+1})$
    $s \leftarrow s'$
  $G_t \leftarrow 0$
  对 $Episode$ 从后向前遍历, $t=\tau-1,\tau-2,...,0$
    取出 $(s_t,r_{t+1})$
    $G_t \leftarrow \gamma G_t+r_{t+1}$
    $G(s_t) \leftarrow G(s_t)+G_t$
    $N(s_t) \leftarrow N(s_t)+1$
$V(S) \leftarrow G(S) / N(S)$
输出:$V(S)$
---
为什么叫做**每次访问法**呢?区别于**首次访问法**,在算法中,在遍历过程中,每次只要遇到 $s_t$ 就要计算它的 $G$ 值而不管它出现在哪里或者出现了多少次。Sutton 认为这两者具有不同的理论基础,首次访问法被研究得很多,具有很强的理论基础;而每次访问法可以过渡到后面要学到的函数近似与资格迹方法,而且从实际效果来看,每次访问法能稍微快一些地收敛,因为样本利用效率更高,相当于采样次数增加了 2~3 倍。
在算法中,这两者的唯一区别就是在对 $Episode$ 从后向前遍历时,是否要检查 $s_t$ 是否曾经出现过。抛开这个区别,算法 10.4 与算法 10.2 的区别是:
- 算法 10.4 要输入策略 $\pi$, 因为要估算 $V_\pi$
- 算法 10.4 要从策略 $\pi$ 选择动作 $a$ 并执行,但是并不会把 $a$ 加入到分幕采样序列中,因为计算 $V_\pi$ 不需要动作信息。
#### 算法实现
【代码位置】Algorithm.Algo_MC_Policy_Evaulation.py
下面我们把代码分成几个小片段来讲解。
函数定义即初始化部分:
```python
# MC 策略评估(预测):每次访问法估算 V_pi
def MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy):
nS = env.observation_space.n # 状态空间
nA = env.action_space.n # 动作空间
Value = np.zeros(nS) # G 的总和
Count = np.zeros(nS) # G 的数量
```
输入:
- env - 环境。
- start_state - 起始状态。
- episodes - 最大循环幕数。
- gamma - 折扣。
- policy - 策略。
初始化部分:
- 获得环境的状态空间,以便初始化数组。
- 获得动作空间,以便在后面选择动作时使用。
- Value 数组保存每个状态上的累积 G 值,除以下面的 Count就会得到平均值 V。
- Count 数组保存每个状态上的累积数量。
多幕循环部分:
```python
for episode in tqdm.trange(episodes): # 多幕循环
Episode = [] # 一幕内的(状态,奖励)序列
s = env.reset() # 重置环境,开始新的一幕采样
done = False
```
这个片段进行多幕循环,为每一幕初始化:
- Episode 保存分幕采样序列。
- s=env.reset() 初始化为起始状态。
- done 用于判断是否结束本幕。
幕内循环部分:
```python
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s]) # 根据策略选择动作
next_s, reward, done, _ = env.step(action) # 与环境交互
Episode.append((s, reward)) # 保存数据序列
s = next_s # 迭代
```
这一片段进行分幕采样,从策略得到动作,执行动作得到 $s',r_{t+1}$,并保存到序列中。当 done 为 True 时表示分幕结束,到达终止状态。
$G$ 值计算部分:
```python
# 从后向前遍历计算 G 值
G = 0 # 不要忘记清零
for t in range(len(Episode)-1, -1, -1):
s, r = Episode[t] # 取出数据
G = gamma * G + r # 计算 G_t
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
```
这一片段进行倒序遍历,计算每个时刻关于 $s$ 的 $G_t$,保存在该 $s$ 的 Value 数组内,计数加 1。
主程序部分:
【代码位置】MC_104_FrozenLake_V_Evaluation.py
```python
if __name__=="__main__":
gamma = 1
episodes = 50000
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=False)
# 随机策略
nA = env.action_space.n
nS = env.observation_space.n
policy = np.ones(shape=(nS, nA)) / nA # 每个状态上的每个动作都有0.25的备选概率
# DP 方法
V_real = get_groud_truth(env, policy, gamma)
# MC 方法
start_state, info = env.reset(seed=5, return_info=True)
V = algoMC.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
env.close()
```
主程序如上,设置参数,初始化环境,生成策略,用 DP 方法得到准确解,再用 MC 方法循环多次做比对,最后关闭环境。
### 10.4.3 误差的平均值 vs. 平均值的误差
很多文献都提到了蒙特卡洛方法是具有高方差低偏差的特点,但是如何证明呢?笔者在这里做了一个实验,来大致说明偏差方差的问题。
假设有一个均值为 3.1 的不知名分布,我们从中采样四次,得到序列 $x=[1,2,4,5]$,则:
- 方差为:$var(x)=[(1-3)^2+(2-3)^2+(4-3)^2+(5-3)^2]/4=2.5$
- 偏差为:$bias(x)=\sqrt{[(1+2+4+5)/4-3.1]^2}=0.1$
方差很大,偏差很小。
根据这个提示,我们可以设计这样一个试验:
1. 循环 2000 幕,根据这 2000 幕的采样数据,计算出 V 值;
2. 做 step 1 十次,就会得到 $V_1,V_2,\cdots,V_{10}$ 共 10 个值;
3. 用动态规划迭代法计算出 $V_*$,作为计算 RMSE 误差的基准;
4. 用 $V_1,V_2,\cdots,V_{10}$ 分别与 $V_*$ 计算 RMSE得到 10 个误差值 $E_1,E_2,\cdots,E_{10}$,然后计算它们的均值 $E_0=(E_1+\cdots+E_{10})/10$,这就是误差的平均值;
5. 计算 $V_0=(V_1+\cdots+V_{10})/10$,得到 $V$ 的平均值,然后计算 $V_0$ 与 $V_*$ 的 RMSE得到 $E'_0$,这就是平均值的误差;
6. 比较 $E_0$ 与 $E'_0$ 的大小。
另外,我们在 2000 幕循环中,每隔 100 次就计算一次上述的 $E_0,E'_0$,这样一共会得到 20 组数据,绘制变化曲线如图 10.4.2。
<center>
<image src="./img/MC-104-V-Error.png">
图 10.4.2 循环 2000 次的误差
左图10 次运行得到的 V 的平均值的误差 $E'_0$右图10 次运行得到的 V 的误差的平均值 $E_0$
</center>
比较左右两图,可以看到左图所示的平均值的误差,起点就比较小,在 0.03,后面的误差可以小到 0.005;右图的误差起点在 0.1 左右,最后的误差是 0.02。所以,$E'_0 < E_0$而且是小一个数量级
- $E_1,\cdots,E_{10}$ 每一个的误差值都很大,说明 $V_1,\cdots,V_{10}$ 的方差较大。
- $V_0$ 是 $V_1,\cdots,V_{10}$ 的均值,虽然它们每一个距离 $V_*$ 都比较远,但是它们的平均值 $V_0$ 距离 $V_*$ 比较近,说明偏差较小。
打印出两组试验的最后一步的数据如下:
```
====================平均值的误差 E'0 ====================
----- 状态价值 V 的平均值 -----
[[0.014 0.012 0.021 0.013]
[0.016 0. 0.041 0. ]
[0.034 0.085 0.145 0. ]
[0. 0.163 0.438 0. ]]
----- 状态价值 V 的平均值的误差 -----
0.0034371081743170516
====================误差的平均值 E0 ====================
----- 最后一个状态价值 V -----
[[0.012 0.008 0.018 0.014]
[0.018 0. 0.033 0. ]
[0.042 0.094 0.146 0. ]
[0. 0.145 0.45 0. ]]
----- 最后一个状态价值 V 的误差 -----
0.01752455053676336
```
有些读者可能会怀疑循环次数不够,造成右图还没有收敛到最佳状态,所以我们把循环次数从 2000 次增加到 10000 次,再看结果如图 10.4.3 所示。
<center>
<image src="./img/MC-104-V-Error2.png">
图 10.4.3 循环 10000 次的误差
左图10 次运行得到的 V 的平均值的误差 $E'_0$右图10 次运行得到的 V 的误差的平均值 $E_0$
</center>
可以看到两者都收敛了,左图的误差可以到 0.002 左右,而右图的误差在 0.015 左右,还是相差了一个数量级。
所以,我们得到的结论是,对于这种高方差低偏差问题,某一次的结果很难让人满意。如何减少方差呢?有兴趣的读者可以自行查询一下以下概念:
- 分层采样法
- 重要采样法
- 相关采样法
- 对偶变量法
当然,最简单的办法,是利用其低偏差的特点,多次试验取平均,往往可以事半功倍。

Просмотреть файл

@ -0,0 +1,212 @@
## 10.5 策略评估(预测)$Q_\pi$
### 10.5.1 策略
在 10.4 小节中,我们学习如何预测 $V_\pi$,而且也提到过,预测 $V_\pi$ 的意义其实不是很大,更有用的还是预测 $Q_\pi$。所以,在本小节中,我们将会使用每次访问法估算指定策略的 $Q_\pi$ 值。
在第八章的射击气球问题中,开始引入了策略,当时只有两个动作:射击红色气球或者蓝色气球,非常简单。第九章的穿越虫洞问题中,我们使用价值迭代直接找到了最优策略,但其实对策略的理解还不够深入。所以在本章中利用冰面行走问题,更深入地研究一下策略问题。
在 10.4 节中,使用了随机策略,即在四个方向上都有 0.25 的概率被选择执行。那么非随机(定向)策略应该如何定义呢?见图 10.5.1。
<center>
<image src="./img/MC-105-3Policy.png">
图 10.5.1 三种策略
(左图:正确策略;中图:随机策略;右图:错误策略)
</center>
在图 10.5.1 中,我们给出了三种策略定义,描述如下:
- 正确策略
因为起点在左上角,终点在右下角,所以给出的策略是各有 0.3 的概率向下方和右侧移动,各有 0.2 的概率向错误方向移动。当然,这里的“错误方向”只是针对这个问题的设置而言,在有些设置中,智能体不得不向左侧或者上方移动来避免掉入冰洞。而“正确方向”也只是在大方向上是正确的,具体到某个状态上,比如 $s_4$,就不能向右侧移动,而是必须向下方甚至向上方移动才是安全的。但是智能体依然会在 $s_4$ 尝试向右侧移动,通过获得惩罚来学习如何找到安全路径。
- 随机策略
向四个方向移动的概率各是 0.25。
- 错误策略
每一步都有相对较大的概率 0.3 向与目标位置相反的方向移动,以较小的概率 0.2 向右下方移动。
用代码来描述的话是这样的:
```python
policy_names = ["正确策略", "随机策略", "错误策略"]
policies = [
#left, down, right, up 四个值相加等于 1
(0.20, 0.30, 0.30, 0.20), # 正确策略
(0.25, 0.25, 0.25, 0.25), # 随机策略
(0.30, 0.20, 0.20, 0.30) # 错误策略
]
```
policies 数组的每一行都是一个策略,由于策略是一种概率表达,所以 4 个值相加要等于 1。
如果定义一个极端的策略,比如 (0.5, 0, 0, 0.5),那么智能体将无法学习到安全路径。
### 10.5.2 算法
接下来要使用每次访问法来预测 $Q$ 值。
#### 算法说明
【算法 10.5】
---
输入:起始状态 $s$,策略 $\pi$,折扣 $\gamma$, 幕数 Episodes
初始化数组:$G(S,A) \leftarrow 0, N(S,A) \leftarrow 0$$S$ 为状态空间,$A$ 为动作空间
多幕 Episodes 循环:
  列表置空 $Episode = [\ ] $ 用于存储序列数据 $(s,a,r)$
  幕内循环直到终止状态:
    从 $s$ 根据策略 $\pi$ 得到动作 $a$
    执行 $a$ 从环境得到 $s',r$ 以及是否终止的标志
    $Episode \Leftarrow (s,a,r)$,相当于是 $(s_t,a_t,r_{t+1})$
    $s \leftarrow s'$
  $G_t \leftarrow 0$
  对 $Episode$ 从后向前遍历, $t=\tau-1,\tau-2,...,0$
    取出 $(s_t,a_t,r_{t+1})$
    $G_t \leftarrow \gamma G_t+r_{t+1}$
    $G(s_t,a_t) \leftarrow G(s_t,a_t)+G_t$
    $N(s_t,a_t) \leftarrow N(s_t,a_t)+1$
$Q(S,A) \leftarrow G(S,A) / N(S,A)$
输出:$Q(S,A)$
---
与算法 10.4 的主要区别:
- 初始化的数组为二维 (S,A),记录每个状态 $s$ 下的每个动作 $a$ 的 G 值。
- Episode 中保存 $s,a,r$,比算法 10.4 多一个动作 $a$。
#### 算法实现
【代码位置】Algo_MC_Policy_Evaulation.py
```python
# MC 策略评估(预测):每次访问法估算 Q_pi
def MC_EveryVisit_Q_Policy(env, start_state, episodes, gamma, policy):
nA = env.action_space.n # 动作空间
nS = env.observation_space.n # 状态空间
Value = np.zeros((nS, nA)) # G 的总和
Count = np.zeros((nS, nA)) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
s = env.reset() # 重置环境,开始新的一幕采样
Episode = [] # 一幕内的(状态,动作,奖励)序列
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s]) # 根据策略选择动作
next_s, reward, done, _ = env.step(action) # 与环境交互得到反馈
Episode.append((s, action, reward)) # 保存 s,a,r
s = next_s
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s,a] += G # 值累加
Count[s,a] += 1 # 数量加 1
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
Q = Value / Count # 求均值
return Q
```
上述代码中与预测 V 函数不同的地方在于 Value 和 Count 数值都是二维的,因为加入了“动作”这一维信息。
### 10.5.3 结果
#### 三种策略的误差比较
以下是主程序中的主要参数设置:
```
gamma = 0.9
episodes = 20000
```
折扣值 0.9,每个策略做 20000 幕采样。
但是我们从 10.4 节学习到,蒙特卡洛法是高偏差的,所以需要把三种策略各运行 10 次(每次都是 20000 幕采样),求其误差的平均值(不是平均值的误差),得到图 10.5.2。
<center>
<image src="./img/MC-105-3-Policy-error.png">
图 10.5.2 三种策略 10 次运行的误差的平均值变化比较
</center>
这个结果可能和读者的想象不太一样:三种策略的误差随着循环次数的增加并没有很大的区别。这说明在不同的策略下,算法都可以评估出正确的动作价值函数,而且收敛速度没有差别。
其它运行指标展示在表 10.5.1 中。
表 10.5.1 三种策略的结果比较
||正确策略|随机策略|错误策略|
|-|-|-|-|
|$\pi(s)$|[0.2, 0.3, 0.3, 0.2]|[0.25, 0.25, 0.25, 0.25]|[0.3, 0.2, 0.2, 0.3]|
|迭代次数(阈值=1e-3)|10|8|6|
|迭代次数(阈值=1e-4)|15|14|12|
|迭代次数(阈值=1e-5)|20|20|20|
|迭代次数(阈值=1e-6)|24|26|28|
|迭代次数(阈值=1e-7)|28|32|35|
|评估时长(秒)|5|7|9|
|平均每幕长度|6.1|7.7|10.4|
数据说明:
- 当阈值为 1e-4 时,用贝尔曼方程的迭代法都是 15 次以内可以收敛,这和图 10.5.2 是相符的,也就是策略不影响收敛速度。而且有趣的是,我们认为是正确的策略,迭代次数反而比错误的策略要多 3 次。
但是当阈值精度提高时,正确策略用的迭代次数就会比错误策略少了。
- 但是用时方面就会有明显差异。单次运行时,正确策略用时最短,错误策略用时最长。
为什么错误策略用时长呢?因为方向性的错误,所以导致智能体需要花费更多的步数达到目的地,正确策略最短,平均用 6.1 步,而错误策略最长,平均用 10.4 步(因为大概率会走反)。
#### 动作价值函数结果
图 10.5.3 显示了三种策略下的 16 个状态上的 4 个动作的价值,笔者特意用红色把 4 个动作的最大值标记出来,便于读者快速阅读。
<center>
<image src="./img/Q-Value.png">
图 10.5.3 三种策略的动作价值函数
(左图:正确策略;中图:随机策略;右图:错误策略)
</center>
读者可以发现,左图中的数值普遍较大,而右图中的数值普遍较小。这是为什么呢?
我们回忆一下第九章中的最优动作价值函数的定义:
$$
\begin{aligned}
q_*(s,a) &\doteq \max_\pi q_\pi(s,a)
\\
&=\max\big[q_{\pi_1}(s,a),q_{\pi_2}(s,a),\cdots,q_{\pi_m}(s,a)\big]
\end{aligned}
\tag{10.5.1}
$$
通过式10.5.1)可以知道,不同的策略 $\pi$ 通过策略评估会得到不同的动作价值函数值 $q_\pi$,值越大的越接近最优动作价值。在图 10.5.3 中,左图是正确的策略(但可能不是最优策略),而右图的策略肯定要比左图的差,所以右图的数值就会小。
但是不管总体值是大是小,在一个策略中的一个状态的四个动作中,总会有一个最大值,比如 $s_0$ 的第二个动作(向下方移动),而三个策略都是在相同位置上具有最大值,所以我们把这个最大值所代表的动作绘制在图 10.5.4 中。
<center>
<image src="./img/Q-MaxValue.png">
图 10.5.4 三种策略的最佳动作方向
</center>
读者会发现有一些状态的动作价值全为 0是因为该状态是终止状态位于 $s_5,s_7,s_{11},s_{12},s_{15}$,包括 4 个冰洞和一个目的地。
读者在这里可能会产生一个问题,为什么看似错误的策略,和一个正确的策略,都能得到同样的正确的结果呢?
我们所说的“错误策略”实际上还有没有错到离谱的程度比如读者可以试试0.5, 0.0, 0.0, 0.5)这个策略,根本得不到正确路径。所以严格地说,上面三个策略并没有正确与错误的区别,只是花费的力气不同,一个具有好的 rule-based基于规则的策略可以帮助智能体更快地学习到收敛状态。
比如,一个经验丰富的驾驶员,都知道有“让车不让道”的说法,意思是遇到前方或侧方的车或行人,第一反应是刹车,而不是躲避,因为躲避有可能带来更大的风险。在无人驾驶问题中,如果有上面这条规则,将会大大提高无人驾驶的安全性。
冰面行走问题是一个简单的问题,我们可以给出所谓的“正确策略”。但是遇到复杂问题时,不能明确什么是正确的,就用随机策略做初始化即可。

Просмотреть файл

@ -0,0 +1,4 @@
动画展示最终过程

Просмотреть файл

@ -0,0 +1,14 @@
### 思考与练习
1. 使用随机投点法计算图 10.4.1 三个函数的定积分值。
2. 10.3.3 节中的冰面行走问题,能用贝尔曼方程的矩阵形式得到解析解吗?为什么?
3. 使用每次访问法估算 10.2 节中的状态价值函数,并于首次访问法比较。
4. 在 10.2 节中关于收敛部分的代码中,增加一些逻辑,来判断一共发生了几次假收敛。
5. 在 10.2 节中,如果运行 10 万次,误差的曲线是什么样子的?看看后续是否很平缓。如果 delta=1e-3, 还能收敛吗?
6. 在 10.5 节中使用一个极端策略比如0.4, 0.1, 0.0, 0.5),看看误差变化情况以及最后是否能得到正确路径。
7. 根据 10.5 节的知识,把 FrozenLake 设置成 8x8然后评估三种策略的 Q 值。
### 参考资料
- Richard S. Sutton 《强化学习》
- https://new.qq.com/omn/20220314/20220314A09IY600.html

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 72 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 57 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 32 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 25 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 8.9 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 11 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 136 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 59 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 53 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 27 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 39 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 45 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 47 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 8.2 KiB

Просмотреть файл

До

Ширина:  |  Высота:  |  Размер: 59 KiB

После

Ширина:  |  Высота:  |  Размер: 59 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 32 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 43 KiB

Просмотреть файл

@ -0,0 +1,61 @@
import numpy as np
import tqdm
# MC 策略评估(预测):每次访问法估算 V_pi
def MC_EveryVisit_V_Policy(env, episodes, gamma, policy):
nS = env.observation_space.n
nA = env.action_space.n
Value = np.zeros(nS) # G 的总和
Count = np.zeros(nS) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
Episode = [] # 一幕内的(状态,奖励)序列
s = env.reset() # 重置环境,开始新的一幕采样
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, _ = env.step(action)
Episode.append((s, reward))
s = next_s
# 从后向前遍历计算 G 值
G = 0
for t in range(len(Episode)-1, -1, -1):
s, r = Episode[t]
G = gamma * G + r
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
V = Value / Count # 求均值
return V
# MC 策略评估(预测):每次访问法估算 Q_pi
def MC_EveryVisit_Q_Policy(env, episodes, gamma, policy):
nA = env.action_space.n # 动作空间
nS = env.observation_space.n # 状态空间
Value = np.zeros((nS, nA)) # G 的总和
Count = np.zeros((nS, nA)) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
# 重置环境,开始新的一幕采样
s, _ = env.reset(return_info=True)
Episode = [] # 一幕内的(状态,动作,奖励)序列
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, _ = env.step(action)
Episode.append((s, action, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s,a] += G # 值累加
Count[s,a] += 1 # 数量加 1
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
Q = Value / Count # 求均值
return Q

Просмотреть файл

@ -0,0 +1,121 @@
import numpy as np
import tqdm
import math
# 来自于第六章
# MC0 - 多次采样并随采样顺序正向计算 G 值,
# 然后获得回报 G 的数学期望,即状态价值函数 v(start_state)
def MC_Sequential_V(dataModel, start_state, episodes, gamma):
G_sum = 0 # 定义最终的返回值G 的平均数
# 循环多幕
for _ in tqdm.trange(episodes):
s = start_state # 把给定的起始状态作为当前状态
G = 0 # 设置本幕的初始值 G=0
t = 0 # 步数计数器
is_end = False
while not is_end:
s_next, reward, is_end = dataModel.step(s)
G += math.pow(gamma, t) * reward
t += 1
s = s_next
# end while
G_sum += G # 先暂时不计算平均值,而是简单地累加
# end for
v = G_sum / episodes # 最后再一次性计算平均值,避免增加计算开销
return v
# MC1-FirstVisit 首次访问法
def MC_FirstVisit_V(dataModel, start_state, episodes, gamma, checkpoint=100, delta=1e-3):
Value = np.zeros(dataModel.nS) # G 的总和
Count = np.zeros(dataModel.nS) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
Episode_State = [] # 一幕内的状态序列
Episode_Reward = [] # 一幕内的奖励序列
s = start_state
is_end = False
while (is_end is False): # 幕内循环
next_s, r, is_end = dataModel.step(s) # 从环境获得下一个状态和奖励
Episode_State.append(s.value)
Episode_Reward.append(r)
s = next_s
assert(len(Episode_State) == len(Episode_Reward))
num_step = len(Episode_State)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s = Episode_State[t]
r = Episode_Reward[t]
G = gamma * G + r
if not (s in Episode_State[0:t]):# 首次访问型
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
return Value / Count # 求均值
# MC - FirstVisit with EarlyStop 首次访问法,早停
def MC_FirstVisit_EarlyStop(dataModel, start_state, episodes, gamma, checkpoint=100, delta=1e-3):
Value = np.zeros(dataModel.nS) # G 的总和
Count = np.zeros(dataModel.nS) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
Episode_State = [] # 一幕内的状态序列
Episode_Reward = [] # 一幕内的奖励序列
s = start_state
is_end = False
while (is_end is False): # 幕内循环
next_s, r, is_end = dataModel.step(s) # 从环境获得下一个状态和奖励
Episode_State.append(s.value)
Episode_Reward.append(r)
s = next_s
assert(len(Episode_State) == len(Episode_Reward))
num_step = len(Episode_State)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s = Episode_State[t]
r = Episode_Reward[t]
G = gamma * G + r
if not (s in Episode_State[0:t]):# 首次访问型
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
if (episode+1)%checkpoint == 0:
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
V = Value / Count # 求均值
if abs(V-V_old).max() < delta: # 判断是否收敛
less_than_delta_count += 1 # 计数器 +1
if (less_than_delta_count == 3): # 计数器到位
break # 停止迭代
else: # 没有持续收敛
less_than_delta_count = 0 # 计数器复位
V_old = V.copy()
return Value / Count # 求均值
# MC2-EveryVisit - 每次访问法
def MC_EveryVisit_V(dataModel, start_state, episodes, gamma):
Value = np.zeros(dataModel.nS) # G 的总和
Count = np.zeros(dataModel.nS) # G 的数量
for _ in tqdm.trange(episodes): # 多幕循环
Episode = [] # 一幕内的(状态,奖励)序列
s = start_state
is_end = False
while (is_end is False): # 幕内循环
next_s, reward, is_end = dataModel.step(s) # 从环境获得下一个状态和奖励
Episode.append((s.value, reward))
s = next_s
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, r = Episode[t]
G = gamma * G + r
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
return Value / Count # 求均值

Просмотреть файл

@ -0,0 +1,57 @@
import numpy as np
import copy
# 式 (9.6.1) 计算 q假设已知 v*(s')
def q_star(p_s_r_d, gamma, V):
q = 0
# 遍历每个转移概率,以计算 q
for p, s_next, reward, done in p_s_r_d:
# math: \sum_{s'} p_{ss'}^a [ r_{ss'}^a + \gamma * v_{\pi}(s')]
q += p * (reward + gamma * V[s_next] * (1-done))
return q
# 式 (9.6.3) 计算 v*
def v_star(s, actions, gamma, V, Q):
list_q = [] # 准备列表记录所有下游的 q*
for action in list(actions.keys()): # actions 是一个字典数据key 是动作
# 遍历每个动作以计算q值进而计算v值
q = q_star(actions[action], gamma, V) # 计算 q*
list_q.append(q) # 加入列表
Q[s,action] = q # 记录下所有的q(s,a)值,不需要再单独计算一次
return max(list_q) # 返回几个q*中的最大值,即 v=max(q)
def calculate_VQ_star(env, gamma, max_iteration):
V = np.zeros(env.observation_space.n) # 初始化 V(s)
Q = np.zeros((env.observation_space.n, env.action_space.n)) # 初始化 Q(s,a)
count = 0
# 迭代
while (count < max_iteration):
V_old = V.copy()
# 遍历所有状态 s
for s in range(env.observation_space.n):
actions = env.P[s]
V[s] = v_star(s, actions, gamma, V, Q)
# 检查收敛性
if abs(V-V_old).max() < 1e-4:
break
count += 1
# end while
print("迭代次数 = ",count)
return V, Q
def get_policy(env, V, gamma):
policy = np.zeros((env.nS, env.nA))
for s in range(env.nS):
actions = env.get_actions(s)
list_q = []
if actions is None:
continue
# 遍历每个策略概率
for action, next_p_s_r in actions:
q_star = 0
for p, s_next, r in next_p_s_r:
q_star += p * (r + gamma * V[s_next])
list_q.append(q_star)
policy[s, np.argmax(list_q)] = 1
return policy

Просмотреть файл

@ -0,0 +1,41 @@
import numpy as np
# 式 (8.4.2) 计算 q_pi
def q_pi(p_s_r_d, gamma, V):
q = 0
# 遍历每个转移概率,以计算 q_pi
for p, s_next, reward, done in p_s_r_d:
# math: \sum_{s'} p_{ss'}^a [ r_{ss'}^a + \gamma * v_{\pi}(s')]
q += p * (reward + gamma * V[s_next])
return q
# 式 (8.4.5) 计算 v_pi
def v_pi(policy, s, actions, gamma, V, Q):
v = 0
for action in list(actions.keys()): # actions 是一个字典数据key 是动作
q = q_pi(actions[action], gamma, V)
# math: \sum_a \pi(a|s) q_\pi (s,a)
v += policy[s][action] * q
# 顺便记录下q(s,a)值,不需要再单独计算一次
Q[s,action] = q
return v
# 迭代法计算 v_pi
def calculate_VQ_pi(env, policy, gamma, iteration, delta=1e-4):
V = np.zeros(env.observation_space.n) # 初始化 V(s)
Q = np.zeros((env.observation_space.n, env.action_space.n)) # 初始化 Q(s,a)
count = 0 # 计数器,用于衡量性能和避免无限循环
# 迭代
while (count < iteration):
V_old = V.copy() # 保存上一次的值以便检查收敛性
# 遍历所有状态 s
for s in range(env.observation_space.n):
actions = env.P[s]
V[s] = v_pi(policy, s, actions, gamma, V, Q)
# 检查收敛性
if abs(V-V_old).max() < delta:
break
count += 1
# end while
print("迭代次数 = ",count)
return V, Q

Просмотреть файл

@ -0,0 +1,105 @@
import numpy as np
import gym
def play_policy(env, policy, render=False):
total_reward = 0.
observation = env.reset(seed=5)
while True:
if render:
env.render() # 此行可显示
action = np.random.choice(env.action_space.n,
p=policy[observation])
observation, reward, done, _ = env.step(action)
total_reward += reward # 统计回合奖励
if done: # 游戏结束
break
return total_reward
def policy_v2q(env,v,s=None,gamma=1.0):
if s is not None: # 计算第v[s]
q = np.zeros(env.action_space.n)
for action in range(env.action_space.n):
for p,next_state,reward,done in env.unwrapped.P[s][action]:
q[action] += p*(reward + gamma*v[next_state]*(1-done))
else:
q = np.zeros((env.observation_space.n,env.action_space.n))
for state in range(env.observation_space.n):
q[state] = policy_v2q(env,v,state,gamma)
return q
# 下面是策略评估
def policy_evluate(env,policy,gamma=1.0,tolerant=1e-4):
v = np.zeros(env.observation_space.n)
vs = np.zeros(env.observation_space.n)
q = np.zeros(env.action_space.n)
while True:
delta = 0
for state in range(env.observation_space.n):
vs[state] = sum(policy_v2q(env,v,state,gamma)*policy[state])
delta = max(delta,vs[state]-v[state])
v[state] = vs[state]
if delta<tolerant :
break
return v
# 下面是策略改进
def policy_improve(env,v,policy,gamma=1.0):
optimal = True
for state in range(env.observation_space.n):
q = policy_v2q(env,v,state,gamma)
a_new = np.argmax(q) # 该步执行贪心策略
if policy[state][a_new] != 1:
policy[state] = 0
policy[state][a_new] = 1
optimal = False
return optimal,policy
# 下面是进行策略迭代
def policy_iterative(env,gamma=1.0,tolerant=1e-6):
# 下面先进行随机策略初始化
policy = np.ones((env.observation_space.n, env.action_space.n)) / env.action_space.n
# 不断进行策略评估与策略改进
while True:
v = policy_evluate(env,policy,gamma,tolerant)
optimal, policy = policy_improve(env,v,policy,gamma)
if optimal == True:
break
return policy,v
if __name__=="__main__":
# 下面尝试用随机策略玩
env = gym.make('FrozenLake-v1', map_name = "4x4", is_slippery=True)
random_policy = \
np.ones((env.observation_space.n, env.action_space.n)) / env.action_space.n
episode_rewards = [play_policy(env, random_policy) for _ in range(100)]
print("随机策略 平均奖励:{:.2f}".format(np.mean(episode_rewards)))
# 下面尝试自己编写的策略评估函数
print('状态价值函数:')
v_random = policy_evluate(env,random_policy)
print(v_random.reshape(4,4))
print('动作函数:')
q_random = policy_v2q(env,v_random,s=None)
print(q_random)
#print('更新前策略是:')
#print(random_policy)
exit(0)
# 下面开始进行策略改进
optimal,policy = policy_improve(env, v_random, random_policy, gamma=1.0)
# 更新前策略是:
print('在一次策略改进之后:')
if optimal == True:
print('找到最优解,策略是:')
print(policy)
else:
print('未找到最优解,策略是:')
print(policy)
# 下面是完整代码
print('*****************下面是完整策略迭代结果********************')
policy,v = policy_iterative(env, gamma=1.0, tolerant=1e-6)
print('找到的最优策略是:{}'.format(policy))
print('找到最优策略是')
print(np.argmax(policy,axis=1).reshape(4,4))

Просмотреть файл

@ -0,0 +1,89 @@
import gym
import Algorithm.Algo_MonteCarlo_MDP as algo
import common.DrawQpi as drawQ
import numpy as np
import common.CommonHelper as helper
ground_truth = np.array([
0.014, 0.012, 0.021, 0.010,
0.016, 0, 0.041, 0,
0.035, 0.088, 0.142, 0,
0, 0.176, 0.439, 0
])
if __name__=="__main__":
gamma = 1
episodes = 50000
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=True)
nA = env.action_space.n
nS = env.observation_space.n
policy = np.ones(shape=(nS, nA)) / nA # 随机策略每个状态上的每个动作都有0.25的备选概率
start_state, info = env.reset(seed=5, return_info=True)
V = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
print(np.reshape(np.round(V,3),(4,4)))
Q = algo.MC_EveryVisit_Q_Policy(env, start_state, episodes, gamma, policy)
print(np.round(Q,3))
drawQ.draw(Q,(4,4))
env.close()
exit(0)
episodes = 50000
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=True)
nA = env.action_space.n
nS = env.observation_space.n
policy = np.ones(shape=(nS, nA)) / nA # 随机策略每个状态上的每个动作都有0.25的备选概率
start_state, info = env.reset(seed=5, return_info=True)
V1 = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
print(np.reshape(np.round(V1,3),(4,4)))
env.close()
print(helper.RMSE(V1, ground_truth))
episodes = 15000
V1 = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
V2 = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
V3 = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
print(np.reshape(V1, (4,4)))
print(helper.RMSE(V1, ground_truth))
print(np.reshape(V2, (4,4)))
print(helper.RMSE(V2, ground_truth))
print(np.reshape(V3, (4,4)))
print(helper.RMSE(V3, ground_truth))
V = (V1 + V2 + V3)/3
print(np.reshape(V, (4,4)))
print(helper.RMSE(ground_truth, V))
env.close()
exit(0)
Q1 = algo.MC_EveryVisit_Q(env, start_state, episodes, gamma)
Q2 = algo.MC_EveryVisit_Q(env, start_state, episodes, gamma)
Q3 = algo.MC_EveryVisit_Q(env, start_state, episodes, gamma)
'''
V1 = algo.MC_EveryVisit_V(env, start_state, episodes, gamma)
V2 = algo.MC_EveryVisit_V(env, start_state, episodes, gamma)
V3 = algo.MC_EveryVisit_V(env, start_state, episodes, gamma)
env.close()
print(np.reshape(V1, (4,4)))
print(np.reshape(V2, (4,4)))
print(np.reshape(V3, (4,4)))
V = (V1 + V2 + V3)/3
print(helper.RMSE(V1, V))
print(helper.RMSE(V2, V))
print(helper.RMSE(V3, V))
'''
Q = (Q1 + Q2 + Q3)/3
drawQ.draw(Q1, (4,4))
drawQ.draw(Q2, (4,4))
drawQ.draw(Q3, (4,4))
print(Q)
drawQ.draw(Q, (4,4))

Просмотреть файл

@ -0,0 +1,89 @@
import gym
import Algorithm.Algo_MonteCarlo_MDP as algo
import common.DrawQpi as drawQ
import numpy as np
import common.CommonHelper as helper
ground_truth = np.array([
0.014, 0.012, 0.021, 0.010,
0.016, 0, 0.041, 0,
0.035, 0.088, 0.142, 0,
0, 0.176, 0.439, 0
])
if __name__=="__main__":
gamma = 1
episodes = 50000
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=False)
nA = env.action_space.n
nS = env.observation_space.n
policy = np.ones(shape=(nS, nA)) / nA # 随机策略每个状态上的每个动作都有0.25的备选概率
start_state, info = env.reset(seed=5, return_info=True)
V = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
print(np.reshape(np.round(V,3),(4,4)))
Q = algo.MC_EveryVisit_Q_Policy(env, start_state, episodes, gamma, policy)
print(np.round(Q,3))
drawQ.draw(Q,(4,4))
env.close()
exit(0)
episodes = 50000
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=True)
nA = env.action_space.n
nS = env.observation_space.n
policy = np.ones(shape=(nS, nA)) / nA # 随机策略每个状态上的每个动作都有0.25的备选概率
start_state, info = env.reset(seed=5, return_info=True)
V1 = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
print(np.reshape(np.round(V1,3),(4,4)))
env.close()
print(helper.RMSE(V1, ground_truth))
episodes = 15000
V1 = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
V2 = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
V3 = algo.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
print(np.reshape(V1, (4,4)))
print(helper.RMSE(V1, ground_truth))
print(np.reshape(V2, (4,4)))
print(helper.RMSE(V2, ground_truth))
print(np.reshape(V3, (4,4)))
print(helper.RMSE(V3, ground_truth))
V = (V1 + V2 + V3)/3
print(np.reshape(V, (4,4)))
print(helper.RMSE(ground_truth, V))
env.close()
exit(0)
Q1 = algo.MC_EveryVisit_Q(env, start_state, episodes, gamma)
Q2 = algo.MC_EveryVisit_Q(env, start_state, episodes, gamma)
Q3 = algo.MC_EveryVisit_Q(env, start_state, episodes, gamma)
'''
V1 = algo.MC_EveryVisit_V(env, start_state, episodes, gamma)
V2 = algo.MC_EveryVisit_V(env, start_state, episodes, gamma)
V3 = algo.MC_EveryVisit_V(env, start_state, episodes, gamma)
env.close()
print(np.reshape(V1, (4,4)))
print(np.reshape(V2, (4,4)))
print(np.reshape(V3, (4,4)))
V = (V1 + V2 + V3)/3
print(helper.RMSE(V1, V))
print(helper.RMSE(V2, V))
print(helper.RMSE(V3, V))
'''
Q = (Q1 + Q2 + Q3)/3
drawQ.draw(Q1, (4,4))
drawQ.draw(Q2, (4,4))
drawQ.draw(Q3, (4,4))
print(Q)
drawQ.draw(Q, (4,4))

Просмотреть файл

@ -1,5 +1,3 @@
import random
from turtle import color
import matplotlib.pyplot as plt
import numpy as np
@ -31,7 +29,8 @@ def put_points(ax, num_total_points):
ax.grid()
if __name__=="__main__":
fig = plt.figure()
np.random.seed(15)
fig = plt.figure(figsize=(8,2.2))
ax = fig.add_subplot(141)
put_points(ax, 100)

Просмотреть файл

@ -3,8 +3,7 @@ import matplotlib.pyplot as plt
import numpy as np
# 随机铺点
# 随机铺点计算 pi 值
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)) # 计算每个点到中心的距离
@ -15,17 +14,18 @@ def put_points(num_total_points):
return pi
if __name__=="__main__":
np.random.seed(15)
pis = []
# 取 1000,1100,1200,...,20000个点做试验
for n in tqdm.trange(1000,20000,100):
pi = 0
# 重复100次求平均
for j in range(100):
pi += put_points(n)
pis.append(pi/100)
pis.append(pi/100) # 求平均
plt.grid()
plt.plot(pis)
plt.plot([0,200],[3.14159265,3.14159265])
plt.plot([0,200],[3.14159265,3.14159265]) # 画基准线
plt.title(str.format("average={0}", np.mean(pis)))
plt.show()
# print(put_points(1000000))

Просмотреть файл

@ -13,13 +13,14 @@ 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
# 求积分 式10.1.1
def integral(f, a, b, n):
v = 0
repeat = 10
for i in range(repeat):
x = np.random.uniform(a, b, size=(n, 1))
x = np.random.uniform(a, b, size=(n, 1)) # 在[a,b]区间上均匀随机取点
y = f(x)
v += np.sum(y) / n * (b-a)
v += np.sum(y) / n * (b-a) # 式(10.1.2)
return v/repeat
def show(ax, f, a, b, n, v):
@ -41,6 +42,7 @@ def show(ax, f, a, b, n, v):
ax.plot(x,y,'o',markersize=1,c='b')
if __name__=="__main__":
np.random.seed(15)
v1 = integral(f1, 0.2, 1, 10000)
print("S1 =",v1)
v2 = integral(f2, 0, 3.1416, 10000)
@ -48,7 +50,7 @@ if __name__=="__main__":
v3 = integral(f3, 0, 1, 10000)
print("S3 =",v3)
fig = plt.figure()
fig = plt.figure(figsize=(9,3))
ax = fig.add_subplot(131)
show(ax, f1, 0.2, 1, 100, v1)
ax = fig.add_subplot(132)

Просмотреть файл

@ -0,0 +1,111 @@
import numpy as np
import MC_102_SafetyDrive_DataModel as env
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus']=False
# MC2 - 反向计算G值记录每个状态的G值首次访问型
def MC_FirstVisit_test(dataModel, start_state, episodes, gamma, checkpoint):
Vs = []
Value = np.zeros(dataModel.nS) # G 的总和
Count = np.zeros(dataModel.nS) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
Ts = [] # 一幕内的状态序列
Tr = [] # 一幕内的奖励序列
s = start_state
is_end = False
while (is_end is False): # 幕内循环
next_s, r, is_end = dataModel.step(s) # 从环境获得下一个状态和奖励
Ts.append(s.value)
Tr.append(r)
s = next_s
assert(len(Ts) == len(Tr))
num_step = len(Ts)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s = Ts[t]
r = Tr[t]
G = gamma * G + r
if not (s in Ts[0:t]):# 首次访问型
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
if (episode+1)%checkpoint == 0:
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
Vs.append(Value / Count) # 求均值
#Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
#Vs.append(Value / Count) # 求均值
return Vs
# MC2 - 反向计算G值记录每个状态的G值每次访问型
def MC_EveryVisit_test(dataModel, start_state, episodes, gamma, checkpoint):
Vs = []
Value = np.zeros(dataModel.nS) # G 的总和
Count = np.zeros(dataModel.nS) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
Ts = [] # 一幕内的状态序列
Tr = [] # 一幕内的奖励序列
s = start_state
is_end = False
while (is_end is False): # 幕内循环
next_s, r, is_end = dataModel.step(s) # 从环境获得下一个状态和奖励
Ts.append(s.value)
Tr.append(r)
s = next_s
num_step = len(Ts)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s = Ts[t]
r = Tr[t]
G = gamma * G + r
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
if (episode+1)%checkpoint == 0:
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
Vs.append(Value / Count) # 求均值
#Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
#Vs.append(Value / Count) # 求均值
return Vs
if __name__=="__main__":
np.random.seed(5)
dataModel = env.DataModel()
episodes = 20000 # 计算 20000 次的试验的均值作为数学期望值
gamma = 1.0
checkpoint = 100
Vs_first = MC_FirstVisit_test(dataModel, dataModel.S.Start, episodes, gamma, checkpoint)
Vs_every = MC_EveryVisit_test(dataModel, dataModel.S.Start, episodes, gamma, checkpoint)
V_groundTruth = env.Matrix(dataModel, 1)
errors_first = []
for v in Vs_first:
error = helper.RMSE(v, V_groundTruth)
errors_first.append(error)
errors_every = []
for v in Vs_every:
error = helper.RMSE(v, V_groundTruth)
errors_every.append(error)
plt.plot(errors_first[10:], label="首次访问法")
plt.plot(errors_every[10:], label="每次访问法")
plt.legend()
plt.grid()
plt.title("首次访问法 vs. 每次访问法")
plt.show()

Просмотреть файл

@ -0,0 +1,58 @@
import numpy as np
import MC_102_SafetyDrive_DataModel as env
import time
import Algorithm.Algo_MC_Value_Evaluation as algo
import common.CommonHelper as helper
def matrix_method(dataModel, gamma):
helper.print_seperator_line(helper.SeperatorLines.long, "矩阵法")
V_truth = env.Matrix(dataModel, gamma)
print(np.round(V_truth, 2))
return V_truth
def sequential_method(dataModel, gamma):
helper.print_seperator_line(helper.SeperatorLines.long, "顺序计算法")
start = time.time()
episodes = 5000 # 计算 5000 次循环的均值作为数学期望值
V = []
for s in dataModel.S: # 遍历每个状态
v = algo.MC_Sequential_V(dataModel, s, episodes, gamma) # 采样计算价值函数
V.append(v) # 保存到字典中
# 打印输出
helper.print_V(dataModel, np.array(V))
end = time.time()
helper.print_seperator_line(helper.SeperatorLines.middle, "耗时")
print("duration =", end-start)
return V
def first_visit_method(dataModel, gamma):
helper.print_seperator_line(helper.SeperatorLines.long, "首次访问法")
start = time.time()
episodes = 20000 # 计算 50000 次的试验的均值作为数学期望值
V = algo.MC_FirstVisit_V(dataModel, dataModel.S.Start, episodes, gamma)
helper.print_V(dataModel, V)
end = time.time()
helper.print_seperator_line(helper.SeperatorLines.middle, "耗时")
print("duration =", end-start)
return V
if __name__=="__main__":
np.random.seed(15)
gamma = 1.0
print("gamma =", gamma)
dataModel = env.DataModel()
V_truth = matrix_method(dataModel, gamma)
V_seq = sequential_method(dataModel, gamma)
helper.print_seperator_line(helper.SeperatorLines.middle, "顺序计算法与矩阵法之间的误差")
print("RMSE =", helper.RMSE(V_seq, V_truth))
V_firstVisit = first_visit_method(dataModel, gamma)
helper.print_seperator_line(helper.SeperatorLines.middle, "首次访问法与矩阵法之间的误差")
print("RMSE =", helper.RMSE(V_firstVisit, V_truth))

Просмотреть файл

@ -0,0 +1,78 @@
import numpy as np
import MC_102_SafetyDrive_DataModel as env
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus']=False
# MC2 - 反向计算G值记录每个状态的G值首次访问型
def MC_FirstVisit_test(dataModel, start_state, episodes, gamma, checkpoint, delta):
Value = np.zeros(dataModel.nS) # G 的总和
Count = np.zeros(dataModel.nS) # G 的数量
V_old = np.zeros(dataModel.nS)
less_than_delta_count = 0
V_history = []
for episode in tqdm.trange(episodes): # 多幕循环
Ts = [] # 一幕内的状态序列
Tr = [] # 一幕内的奖励序列
s = start_state
is_end = False
while (is_end is False): # 幕内循环
next_s, r, is_end = dataModel.step(s) # 从环境获得下一个状态和奖励
Ts.append(s.value)
Tr.append(r)
s = next_s
assert(len(Ts) == len(Tr))
num_step = len(Ts)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s = Ts[t]
r = Tr[t]
G = gamma * G + r
if not (s in Ts[0:t]):# 首次访问型
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
if (episode+1)%checkpoint == 0:
Count[Count==0] = 1 # 把分母为0的填成1主要是终止状态
V = Value / Count # 求均值
V_history.append(V)
#print(np.round(V,3))
if abs(V-V_old).max() < delta: # 判定收敛
less_than_delta_count += 1 # 计数器 +1
if (less_than_delta_count == 3): # 计数器到位
break # 停止迭代
else: # 没有持续收敛
less_than_delta_count = 0 # 计数器复位
V_old = V.copy()
print("循环幕数 =",episode+1)
return V_history
if __name__=="__main__":
np.random.seed(5)
dataModel = env.DataModel()
episodes = 100000 # 计算 20000 次的试验的均值作为数学期望值
gamma = 1.0
checkpoint = 200
delta = 1e-2
V_history = MC_FirstVisit_test(dataModel, dataModel.S.Start, episodes, gamma, checkpoint, delta)
V_truth = env.Matrix(dataModel, gamma)
Errors = []
for V in V_history:
error = helper.RMSE(V, V_truth)
Errors.append(error)
plt.plot(Errors)
plt.grid()
plt.xlabel("循环次数(x200)")
plt.ylabel("RMSE 误差")
plt.show()
print("状态值 =", np.round(V,3))
print("误差 =", error)

Просмотреть файл

@ -1,7 +1,8 @@
from asyncio import create_subprocess_shell
import numpy as np
from enum import Enum
import common.CommonHelper as helper
# 状态
class States(Enum):
@ -11,7 +12,7 @@ class States(Enum):
DownSpeed = 3 # 闹市减速
ExceedSpeed = 4 # 超速行驶
RedLight = 5 # 路口闯灯
LowSpeed = 6 # 小区减速
LowSpeed = 6 # 小区慢行
MobilePhone = 7 # 拨打手机
Crash = 8 # 发生事故
Goal = 9 # 安全抵达
@ -19,18 +20,18 @@ class States(Enum):
# 状态转移概率
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]
[ #S N P D E R L M C G E
[0.0, 0.9, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], # S
[0.0, 0.0, 0.2, 0.1, 0.1, 0.1, 0.3, 0.1, 0.0, 0.1, 0.0], # N
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], # P
[0.0, 0.7, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], # D
[0.0, 0.2, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.5, 0.0, 0.0], # E
[0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0], # R
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0], # L
[0.0, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0], # M
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0], # C
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0], # G
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0] # E
]
)
@ -43,7 +44,7 @@ class DataModel(object):
self.P = P # 状态转移矩阵
self.R = R # 奖励
self.S = States # 状态集
self.num_states = len(self.S) # 状态数量
self.nS = len(self.S) # 状态数量
self.end_states = [self.S.End] # 终止状态集
self.V_ground_truth = Matrix(self, 1)
@ -57,18 +58,14 @@ class DataModel(object):
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 step(self, s):
next_s = np.random.choice(self.S, p=self.P[s.value])
return next_s, self.get_reward(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.nS) * (1+1e-7)
#I = np.eye(dataModel.num_states)
tmp1 = I - gamma * dataModel.P
tmp2 = np.linalg.inv(tmp1)
@ -77,5 +74,5 @@ def Matrix(dataModel, gamma):
if __name__=="__main__":
dataModel = DataModel()
v = Matrix(dataModel, 1.0)
print(np.around(v,2))
V = Matrix(dataModel, 1.0)
helper.print_V(dataModel, V)

Просмотреть файл

@ -0,0 +1,35 @@
# 用 DP 迭代方法计算随机策略下的 FrozenLake 问题的 V_pi, Q_pi
import numpy as np
import gym
import Algorithm.Algo_PolicyValueFunction as algoPI
import Algorithm.Algo_OptimalValueFunction as algoStar
import common.DrawQpi as drawQ
import common.CommonHelper as helper
if __name__=="__main__":
slips = [False, True]
gamma = 1
iteration = 1000
for slip in slips:
helper.print_seperator_line(helper.SeperatorLines.long, str.format("is slippery={0}", slip))
env = gym.make('FrozenLake-v1', map_name = "4x4", is_slippery=slip)
# print(env.P)
# 随机策略
policy = np.ones((env.observation_space.n, env.action_space.n)) / env.action_space.n
V_pi, Q_pi = algoPI.calculate_VQ_pi(env, policy, gamma, iteration)
helper.print_seperator_line(helper.SeperatorLines.middle, "随机策略下的状态价值函数 V")
print(np.reshape(np.round(V_pi,3), (4,4)))
helper.print_seperator_line(helper.SeperatorLines.middle, "随机策略下的动作价值函数 Q")
print(np.round(Q_pi,3))
drawQ.draw(Q_pi, (4,4))
V_star, Q_star = algoStar.calculate_VQ_star(env, gamma, iteration)
helper.print_seperator_line(helper.SeperatorLines.middle, "最优策略下的状态价值函数 V")
print(np.reshape(np.round(V_star,3), (4,4)))
helper.print_seperator_line(helper.SeperatorLines.middle, "最优策略下的动作价值函数 Q")
Q_star = np.round(Q_star,3)
print(Q_star)
drawQ.draw(Q_star, (4,4))
env.close()

Просмотреть файл

@ -0,0 +1,116 @@
import numpy as np
import gym
import Algorithm.Algo_MC_Policy_Evaulation as algoMC
import Algorithm.Algo_PolicyValueFunction as algoDP
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus']=False
def get_groud_truth(env, policy, gamma):
iteration = 100
V, _ = algoDP.calculate_VQ_pi(env, policy, gamma, iteration)
return V
# MC 策略评估(预测):每次访问法估算 V_pi
def MC_EveryVisit_V_Policy_test(env, episodes, gamma, policy, checkpoint=1000):
nS = env.observation_space.n
nA = env.action_space.n
Value = np.zeros(nS) # G 的总和
Count = np.zeros(nS) # G 的数量
V_old = np.zeros(nS)
V_history = [] # 测试用
for episode in tqdm.trange(episodes): # 多幕循环
Episode = [] # 保存一幕内的(状态,奖励)序列
s, _ = env.reset(return_info=True)# 重置环境,开始新的一幕采样
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, info = env.step(action)
Episode.append((s, reward))
s = next_s
# 从后向前遍历计算 G 值
G = 0
for t in range(len(Episode)-1, -1, -1):
s, r = Episode[t]
G = gamma * G + r
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
# 检查是否收敛
if (episode + 1)%checkpoint == 0:
Count[Count==0] = 1 # 把分母为0的填成1主要是对终止状态
V = Value / Count
V_history.append(V)
return V_history # 返回历史数据用于评测
if __name__=="__main__":
gamma = 1
episodes = 10000
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=False)
# 随机策略
policy = helper.create_policy(
env.observation_space.n, env.action_space.n, (0.25, 0.25, 0.25, 0.25))
# DP
V_truth = get_groud_truth(env, policy, gamma)
# MC
start_state, info = env.reset(seed=5, return_info=True)
# V = algoMC.MC_EveryVisit_V_Policy(env, start_state, episodes, gamma, policy)
checkpoint = 100
runs = 10
count = int(episodes/checkpoint)
V_sum = [] # 准备数据计算平均值的误差
Errors_history = [] # 准备数据计算误差的平均值
for i in range(runs): # 运行 n 次的平均值
V_sum.append([])
Errors_history.append([]) # 增加一维空列表,用于存储历史误差数据
V_history = MC_EveryVisit_V_Policy_test(env, episodes, gamma, policy, checkpoint=checkpoint)
for j in range(len(V_history)):
V = V_history[j]
error = helper.RMSE(V, V_truth)
Errors_history[i].append(error)
V_sum[i].append(V) # 同时计算 V 的平均值
env.close()
# n 次 run 的在每个checkpoint上的平均值的误差
V_average = np.array(V_sum).mean(axis=0)
rmses = np.zeros(count)
for i in range(count):
rmses[i] = helper.RMSE(V_average[i], V_truth)
# 显示平均值的误差
fig = plt.figure(figsize=(10,4))
ax1 = fig.add_subplot(121)
ax1.plot(rmses)
ax1.set_xlabel(u"循环次数(x1000)")
ax1.set_ylabel(u"RMSE 误差")
ax1.set_title(u'平均值的误差')
ax1.grid()
helper.print_seperator_line(helper.SeperatorLines.long, "平均值的误差")
helper.print_seperator_line(helper.SeperatorLines.short, "状态价值 V 的平均值")
print(np.around(V_average[-1], 3).reshape(4,4))
helper.print_seperator_line(helper.SeperatorLines.short, "状态价值 V 的平均值的误差")
print(rmses[-1])
# n 次 run 的在每个checkpoint上的误差的平均值
Errors = np.array(Errors_history).mean(axis=0)
helper.print_seperator_line(helper.SeperatorLines.long, "误差的平均值")
helper.print_seperator_line(helper.SeperatorLines.short, "最后一个状态价值 V")
print(np.around(V,3).reshape(4,4))
helper.print_seperator_line(helper.SeperatorLines.short, "最后一个状态价值 V 的误差")
print(Errors[-1])
# 显示误差的平均值
ax2 = fig.add_subplot(122)
ax2.plot(Errors)
ax2.set_title(u'误差的平均值')
ax2.set_xlabel(u'循环次数(x1000)')
ax2.set_ylabel(u'RMSE 误差')
ax2.grid()
plt.show()

Просмотреть файл

@ -0,0 +1,111 @@
import numpy as np
import gym
import Algorithm.Algo_MC_Policy_Evaulation as algoMC
import Algorithm.Algo_PolicyValueFunction as algoDP
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
def test_n_times(env, episodes, gamma, policy, checkpoint=1000, n=10):
Errors = []
for i in range(n):
Q_history, T_len = MC_EveryVisit_Q_Policy_test(env, episodes, gamma, policy, checkpoint=checkpoint)
errors = []
for Q in Q_history:
error = helper.RMSE(Q, Q_real)
errors.append(error)
Errors.append(errors)
E_array = np.array(Errors)
E_mean = np.mean(E_array, axis=0) # 求 n 次的误差平均值
return E_mean, T_len, Q, error
# MC 策略评估(预测):每次访问法估算 Q_pi
def MC_EveryVisit_Q_Policy_test(env, episodes, gamma, policy, checkpoint=1000):
nA = env.action_space.n
nS = env.observation_space.n
Value = np.zeros((nS, nA)) # G 的总和
Count = np.zeros((nS, nA)) # G 的数量
Q_history = []
T_len = 0
for episode in tqdm.trange(episodes): # 多幕循环
# 重置环境,开始新的一幕采样
s, _ = env.reset(return_info=True)
Episode = [] # 一幕内的(状态,奖励)序列
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, _ = env.step(action)
Episode.append((s, action, reward))
s = next_s
T_len += len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(len(Episode)-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s,a] += G # 值累加
Count[s,a] += 1 # 数量加 1
if (episode + 1)%checkpoint == 0:
Count[Count==0] = 1 # 把分母为0的填成1主要是对终止状态
Q = Value / Count
Q_history.append(Q)
return Q_history, T_len/episodes
def get_groud_truth(env, policy, gamma):
iteration = 100
_, Q = algoDP.calculate_VQ_pi(env, policy, gamma, iteration, delta=1e-7)
return Q
if __name__=="__main__":
gamma = 0.9
episodes = 20000
policy_names = ["正确策略", "随机策略","错误策略"]
policies = [
# left, down, right, up
(0.2, 0.3, 0.3, 0.2),
(0.25, 0.25, 0.25, 0.25),
(0.3, 0.2, 0.2, 0.3)
]
end_states = [5, 7, 11, 12, 15]
np.set_printoptions(suppress=True)
for i, policy_data in enumerate(policies):
helper.print_seperator_line(helper.SeperatorLines.long, policy_names[i])
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=False)
nA = env.action_space.n
nS = env.observation_space.n
policy = helper.create_policy(nS, nA, policy_data)
helper.print_seperator_line(helper.SeperatorLines.short, "策略")
print(policy)
Q_real = get_groud_truth(env, policy, gamma)
start_state, info = env.reset(seed=5, return_info=True)
Errors, T_len, Q, error = test_n_times(env, episodes, gamma, policy, checkpoint=1000, n=1) # n=10
helper.print_seperator_line(helper.SeperatorLines.short, "动作价值函数")
print(np.round(Q,3))
helper.print_seperator_line(helper.SeperatorLines.short, "误差")
print("RMSE =", error)
helper.print_seperator_line(helper.SeperatorLines.short, "平均每幕长度")
print("Len =", T_len)
plt.plot(Errors, label=policy_names[i])
env.close()
policy = helper.extract_policy_from_Q(Q, end_states)
helper.print_seperator_line(helper.SeperatorLines.short, "抽取出来的策略")
print(policy)
drawQ.draw(policy,(4,4))
plt.title(u'三种策略的误差收敛情况比较')
plt.xlabel(u'循环次数(x1000)')
plt.ylabel(u'误差 RMSE')
plt.legend()
plt.grid()
plt.show()

Просмотреть файл

@ -0,0 +1,110 @@
import numpy as np
import gym
import Algorithm.Algo_MC_Policy_Evaulation as algoMC
import Algorithm.Algo_PolicyValueFunction as algoDP
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
def test_n_times(env, episodes, gamma, policy, checkpoint=1000, n=10):
Errors = []
for i in range(n):
Q_history, T_len = MC_EveryVisit_Q_Policy_test(env, episodes, gamma, policy, checkpoint=checkpoint)
errors = []
for Q in Q_history:
error = helper.RMSE(Q, Q_real)
errors.append(error)
Errors.append(errors)
E_array = np.array(Errors)
E_mean = np.mean(E_array, axis=0) # 求 n 次的误差平均值
return E_mean, T_len, Q, error
# MC 策略评估(预测):每次访问法估算 Q_pi
def MC_EveryVisit_Q_Policy_test(env, episodes, gamma, policy, checkpoint=1000):
nA = env.action_space.n
nS = env.observation_space.n
Value = np.zeros((nS, nA)) # G 的总和
Count = np.zeros((nS, nA)) # G 的数量
Q_history = []
T_len = 0
for episode in tqdm.trange(episodes): # 多幕循环
# 重置环境,开始新的一幕采样
s = env.reset()
Episode = [] # 一幕内的(状态,奖励)序列
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, _ = env.step(action)
Episode.append((s, action, reward))
s = next_s
T_len += len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(len(Episode)-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s,a] += G # 值累加
Count[s,a] += 1 # 数量加 1
if (episode + 1)%checkpoint == 0:
Count[Count==0] = 1 # 把分母为0的填成1主要是对终止状态
Q = Value / Count
Q_history.append(Q)
return Q_history, T_len/episodes
def get_groud_truth(env, policy, gamma):
iteration = 100
_, Q = algoDP.calculate_VQ_pi(env, policy, gamma, iteration)
return Q
if __name__=="__main__":
gamma = 0.9
episodes = 20000
policy_names = ["正确策略", "随机策略","错误策略"]
policies = [
# left, down, right, up
(0.2, 0.3, 0.3, 0.2),
(0.25, 0.25, 0.25, 0.25),
(0.3, 0.2, 0.2, 0.3)
]
end_states = [19, 29, 35, 41, 42, 46, 49, 52, 54, 59, 63]
np.set_printoptions(suppress=True)
for i, policy_data in enumerate(policies):
helper.print_seperator_line(helper.SeperatorLines.long, policy_names[i])
env = gym.make("FrozenLake-v1", desc=None, map_name = "8x8", is_slippery=False)
nA = env.action_space.n
nS = env.observation_space.n
policy = helper.create_policy(nS, nA, policy_data)
helper.print_seperator_line(helper.SeperatorLines.short, "策略")
print(policy)
Q_real = get_groud_truth(env, policy, gamma)
start_state, info = env.reset(seed=5, return_info=True)
Errors, T_len, Q, error = test_n_times(env, episodes, gamma, policy, checkpoint=1000, n=1) # n=10
helper.print_seperator_line(helper.SeperatorLines.short, "动作价值函数")
print(np.round(Q,3))
helper.print_seperator_line(helper.SeperatorLines.short, "误差")
print("RMSE =", error)
helper.print_seperator_line(helper.SeperatorLines.short, "平均每幕长度")
print("Len =", T_len)
plt.plot(Errors, label=policy_names[i])
env.close()
policy = helper.extract_policy_from_Q(Q, end_states)
helper.print_seperator_line(helper.SeperatorLines.short, "抽取出来的策略")
print(policy)
drawQ.draw(policy,(8,8))
plt.title(u'三种策略的误差收敛情况比较')
plt.xlabel(u'循环次数(x1000)')
plt.ylabel(u'误差 RMSE')
plt.legend()
plt.grid()
plt.show()

Просмотреть файл

@ -0,0 +1,83 @@
from enum import Enum
import numpy as np
class SeperatorLines(Enum):
empty = 0 # 打印空行
short = 1 # 打印10个'-'
middle = 2 # 打印30个'-'
long = 3 # 打印40个'='
def print_seperator_line(style: SeperatorLines, info=None):
if style == SeperatorLines.empty:
print("")
elif style == SeperatorLines.short:
if (info is None):
print("-"*10)
else:
print("----- " + info + " -----")
elif style == SeperatorLines.middle:
if (info is None):
print("-"*30)
else:
print("-"*15 + info + "-"*15)
elif style == SeperatorLines.long:
if (info is None):
print("="*40)
else:
print("="*20 + info + "="*20)
def print_V(dataModel, V):
vv = np.around(V,2)
print("状态价值函数计算结果(数组) :", vv)
for s in dataModel.S:
print(str.format("{0}:\t{1}", s.name, vv[s.value]))
def RMSE(x, y):
err = np.sqrt(np.sum(np.square(x - y))/y.shape[0])
return err
def test_policy(env, policy, episodes=100):
R = 0
for i in range(episodes):
s = env.reset()
done = False
while (done is False): # 幕内循环
action = np.argmax(policy[s])
# action = np.random.choice(env.action_space.n, p=policy[s])
next_s, reward, done, info = env.step(action)
R += reward
if done == True and reward == 0:
print(s, action, next_s)
s = next_s
return R
# 根据输入的4个概率值创建策略
def create_policy(nS, nA, args):
assert(nA == len(args))
left = args[0]
down = args[1]
right = args[2]
up = args[3]
assert(left+down+right+up==1)
policy = np.zeros((nS, nA))
policy[:, 0] = left
policy[:, 1] = down
policy[:, 2] = right
policy[:, 3] = up
return policy
# 从Q函数表格中抽取策略
def extract_policy_from_Q(Q, end_states):
policy = np.zeros_like(Q)
for s in range(Q.shape[0]):
if s not in end_states:
max_v = np.max(Q[s])
for a in range(Q[s].shape[0]):
if Q[s,a] == max_v:
policy[s, a] = 1
return policy

Просмотреть файл

@ -0,0 +1,74 @@
import numpy as np
LEFT, DOWN, RIGHT, UP = 0, 1, 2, 3
LEFT_ARROW = 0x25c4
UP_ARROW = 0x25b2
RIGHT_ARROW = 0x25ba
DOWN_ARROW = 0x25bc
EMPTY_SPACE = 0x0020
CENTER_CROSS = 0x253c
SEP_LINE = 0x2500
class GridCell(object):
def __init__(self, q):
self.space = np.zeros((3,5), dtype=int)
self.space.fill(EMPTY_SPACE) # 空格
self.space[1,2] = CENTER_CROSS # 中心
self.q = np.round(q, 4)
# policy 是一个1x4的数组或矩阵,如[0.1,0.1,0.4,0.4]
# 这种情况要在向上和向右的地方画两个箭头
# 01234
# 0 ^
# 1 o >
# 2
if np.sum(q) != 0:
best_actions = np.argwhere(self.q == np.max(self.q))
pos = best_actions.flatten().tolist()
for action in pos:
if action == LEFT: # left
self.space[1,0] = LEFT_ARROW
self.space[1,1] = SEP_LINE
elif action == UP: # up
self.space[0,2] = UP_ARROW
elif action == RIGHT: # right
self.space[1,3] = SEP_LINE
self.space[1,4] = RIGHT_ARROW
elif action == DOWN: # down
self.space[2,2] = DOWN_ARROW
class Grid(object):
def __init__(self, Q, shape):
self.array = np.zeros((shape[0]*3, shape[1]*5), dtype=int)
for i in range(len(Q)):
row = (int)(i / shape[0])
col = (int)(i % shape[0])
q = Q[i]
cell = GridCell(q)
self.array[row*3:row*3+3, col*5:col*5+5] = cell.space
def draw(Q, shape):
grid = Grid(Q, shape)
for j, rows in enumerate(grid.array):
if (j % 3 == 0):
print("+-----"*shape[1], end="")
print("+")
print("|", end="")
for i,col in enumerate(rows):
print(chr(col), end="")
if ((i+1) % 5 == 0):
print("|", end="")
print()
print("+-----"*shape[1])
if __name__=="__main__":
Q = np.array([
[0.0155, 0.0164, 0.012, 0.015],
[0.0, 0.0, 0.00, 0.00],
[5.07, 3.06, 7.86 , 2.07],
[8.73, 8.73, 8.73 , 8.73]
])
draw(Q, (2,2))

Просмотреть файл

@ -1,229 +0,0 @@
## 蒙特卡洛方法
本节中将要学习的蒙特卡洛方法,是强化学习的一种重要手段,在后面的学习中,要更充分地讨论蒙特卡洛预测和蒙特卡洛控制,都是以本节的内容为理论基础的。
其实在“三门问题”和“学生学习问题”中,已经使用了蒙特卡洛的思想来解决问题了,接下来我们系统地了解一下这种方法。
### 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">
图 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值也不太容易。
通常蒙特卡洛方法可以粗略地分成两类:
- 所求解的问题本身具有内在的随机性,借助计算机的运算能力可以直接**模拟随机的过程**
例如在核物理研究中,分析中子在反应堆中的传输过程。中子与原子核作用受到量子力学规律的制约,人们只能知道它们相互作用发生的概率,却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向。
科学家依据其概率进行随机抽样得到裂变位置、速度和方向,这样模拟大量中子的行为后,经过统计就能获得中子传输的范围,作为反应堆设计的依据。
- 所求解问题可以转化为某种**随机分布的特征数**
比如随机事件出现的概率,或者随机变量的期望值。通过随机抽样的方法,以随机事件出现的频率估计其概率,或者以抽样的数字特征估算随机变量的数字特征,并将其作为问题的解这种方法多用于求解复杂的多维积分问题。
我们在上面举的例子是一个一维的积分问题,积分结果是曲线下的面积;当然可以再增加一维随机采样,变成二维积分问题,结果为曲面下的体积,比一维积分更有实际作用。
蒙特卡洛在强化学习中应用的核心主要包含以下几点:
- MC方法是直接从经验轨迹当中直接进行学习。
- MC方法是一种model-free方法即没有MDP的转移概率以及奖励的先验知识。
- MC方法从完整的经验轨迹中学习不使用bootstrapping方法。
- MC方法简单得使用这个思想价值=平均回报。
- 缺点在于只能应用于一定有终结点的按幕分的MDP过程。
### 参考资料
https://new.qq.com/omn/20220314/20220314A09IY600.html

Просмотреть файл

@ -1,145 +0,0 @@
## 每次访问型 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">
图 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=\tau-1,\tau-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">
图 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

Просмотреть файл

@ -1,7 +0,0 @@
### 批量更新算法
$$
V(s) = V(s) + \alpha (G_t - V)
$$

Двоичный файл не отображается.

До

Ширина:  |  Высота:  |  Размер: 72 KiB

Двоичный файл не отображается.

До

Ширина:  |  Высота:  |  Размер: 53 KiB

Двоичный файл не отображается.

До

Ширина:  |  Высота:  |  Размер: 9.4 KiB

Двоичный файл не отображается.

До

Ширина:  |  Высота:  |  Размер: 8.0 KiB

Двоичный файл не отображается.

До

Ширина:  |  Высота:  |  Размер: 49 KiB

Двоичный файл не отображается.

До

Ширина:  |  Высота:  |  Размер: 38 KiB

Двоичный файл не отображается.

До

Ширина:  |  Высота:  |  Размер: 22 KiB

Просмотреть файл

@ -1,82 +0,0 @@
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,146 @@
## 11.1 二十一点问题
二十一点是一个十分具有趣味性的牌类游戏最早出现在16世纪起源于法国法语称为ving-et-un20 和 1。后传入英国并广泛流传。如果玩家拿到黑心 A 和黑心 J 会给予额外的奖励英文名称叫黑杰克Blackjack
<center>
<img src="./img/BlackJack.png">
图 11.1.2 二十一点游戏
</center>
该游戏由 2 到 6 人进行,一名庄家,其余为玩家。
使用除大小王之外的 52 张扑克牌,一般为多副牌(比如四副)混在一起使用,以避免依赖一些简单的概率计算来获胜。比如,如果只使用一副牌,则一共有四张 A如果有四个玩家都有一张 A那么庄家肯定没有 A。如果是四副牌的话就不用考虑这些问题了相当于牌池足够大。
玩家的目标是使手中牌的点数之和不超过 21 点且尽量大。
游戏过程中会给庄家和玩家每人发两张牌,一张为明牌,一张为暗牌。接下来会有几种情况:
- 如果玩家手中的两张牌之点数和为 21称为天和除非庄家也是天和算平局否则玩家直接获胜。
- 玩家可以根据手中牌点数和的大小与庄家手中牌点数来要牌 (hit) 并决定何时停牌 (stick) 。当选择继续要牌后若总和大于21点则算自爆bust游戏失败无论后续庄家如何。
- 全部停牌后,庄家翻开扣着的牌,如果点数之和小于 17 点,则继续要牌,直到所有点数之和是 17 点或大于 17 点后,和玩家进行比较,谁的点数更靠近 21谁获胜如果庄家自爆玩家获胜若两方点数相同则为平局。
具体点数计算规则如下:
1. 2 到 10 的点数就是其牌面的数字;
2. J,Q,K 三种牌均记为 10 点;
3. AAce牌可以当作 1 点,也可以当作 11 点11 点时称为“有可用的 A”Usable Ace
#### 动作空间
我们要研究的是玩家的策略,所以只有两个动作:
- 要牌hit= 1
- 停牌stick= 0
#### 状态空间
表 11.1.1 玩家的状态
|庄家手中的明牌 $\to$<br>玩家手中的牌$\downarrow$|1|2|...|10|
|-|-|-|-|-|
|12||||
|13||||
|...||||
|20||||
|21||||
玩家的状态由两个因素决定:
1. 庄家手中的牌(一张明牌);
2. 玩家自己手中的所有牌。
因为当玩家手中的牌不到 12 点时,再来任何一张牌都后,总数都不会超过 21 点,所以我们只考虑 12 点到 21 点(包含天和的情况)之间的情况。
所以,所有的状态在表 11.1.1 中所示,一共有 10x10=100 种状态。
另外,玩家有 A 的话,那将会有很大的变数,所以有 A 无 A 算是第三个因素,那么整个状态空间是 2x10x10=200。
#### 分幕与奖励
- 玩家赢R=1
- 玩家输R=-1
- 平局R=0。
每一局看作一幕,中间过程没有任何奖励。
在这个强化学习问题中,拿到黑桃 A 和 J 不再有奖励,因为它与强化学习无关,只与运气有关。
#### 强化学习的目标
另外,在实际的游戏中,庄家面对多个玩家,假设一个玩家自爆,不影响其它玩家继续,庄家会得到针对该玩家的奖励,即使庄家后面输给其它玩家或自爆。
在本问题中,是站在某个玩家的角度来考虑问题,不管其它玩家的情况如何。因为庄家的规则是固定的,玩家的策略相对灵活,而不同的玩家有不同的策略,因人而异。我们的目的是找到最优的玩家策略。
所以,就假设只有庄家和一个玩家。
#### 环境测试
```python
import gym
import numpy as np
env = gym.make('Blackjack-v1', sab=True) # 使用 Sutton 书中的设置
for i in range(10): # 循环 10 幕
s = env.reset() # 重置环境,开始一幕
Episode = [] # 保存记录
while True:
action = np.random.choice(env.action_space.n) # 随机选择动作做测试
next_s, r, done, _ = env.step(action) # 与环境交互
Episode.append((s, action, r, next_s)) # 保存 (S,A,R,S')
if done: # 本幕结束
break
s = next_s # 迭代
print(Episode)
env.close()
```
其中的参数Argumentssab=True 的含义是使用 Sutton and Barto 书中的设置。
得到的部分典型结果如下:
```
# s, a, r, s'
[((20, 7, False), 1, -1.0, (29, 7, False))]
......
[((19, 2, True), 0, 1.0, (19, 2, True))]
......
[((17, 7, True), 1, 0.0, (18, 7, True)), ((18, 7, True), 1, 0.0, (20, 7, True)),
((20, 7, True), 1, 0.0, (21, 7, True)), ((21, 7, True), 0, 1.0, (21, 7, True))]
......
[((16, 9, False), 1, 0.0, (18, 9, False)), ((18, 9, False), 0, -1.0, (18, 9, False))]
```
- 第一个结果
20初始玩家拿到两张牌 20 点;
7庄家明牌 7 点;
False玩家没有 A
1选择要牌。
-1.0:玩家输了,奖励 -1.0。
原因在 s' 中,三张牌 29 点爆了bust
- 第二个结果
初始玩家拿到两张牌 19 点,庄家明牌 2 点,玩家有 A选择停牌。
赢了,奖励 1.0。原因是庄家可能爆了,或者点数少。
- 第三个结果
初始玩家拿到两张牌 17 点,庄家明牌 7 点,玩家有 A选择要牌
又拿到一张 A18 点,选择要牌;
拿到一张 220 点,选择要牌;
又拿到一张 A太幸运了21 点,选择停牌。
赢了,奖励 1.0。原因是庄家点数小。
- 第四个结果
初始玩家拿到两张牌 16 点,庄家明牌 9 点,玩家没有 A选择要牌
拿到一张 218 点,选择停牌。
输了,奖励 -1.0。原因是庄家点数大。
在结果中并没有说明输或赢的原因,是笔者根据基本逻辑猜测的,对于强化学习问题来说并不重要。

Просмотреть файл

@ -0,0 +1,226 @@
## 11.2 蒙特卡洛控制
什么是蒙特卡洛控制?
策略迭代
上一节中,通过对三种策略的评估比较得知:相对正确的策略,其动作价值函数值也会相对较大。那么如何找到最优策略呢?
我们先复习一下第九章中学习的价值迭代的工作方法,如图 11.1.1 所示。
<center>
<img src='./img/ValueIteration.png'>
图 11.1.1 价值迭代方法示意图
</center>
价值迭代,是以计算动作价值函数 $Q_{k+1}(s,a)=P[R+\gamma V_k(s')]$ 为手段,在每次迭代中都取 $V_{k+1}(s)=\max_a Q_{k+1}(s,a)$,这样一步步地逼近最优价值函数。
在本章中,由于没有环境信息(缺乏 $P,R$ 的直接定义),所以无法使用价值迭代法,只能尝试其它方法。而在第十章中曾经提到过:预测 $Q$ 函数比预测 $V$ 函数更有用,因为从 $Q$ 函数的表格型结果中可以抽取出策略 $\pi$ 来。由此,我们看到了解决问题的思路:
1. 首先,给定任意策略 $\pi$,使用蒙特卡洛方法来预测 $Q$,这在第十章中已经学习过了;
2. 其次,有了 $Q$ 值后,按照价值迭代的办法,从中直接取最大值所对应的动作,由此可以确定新的策略;
3. 根据新的策略,回到第 1 步。
何时结束上述循环过程呢?可以有两个方法:
1. 指定循环次数;
2. 比较旧策略和新策略,如果相等,说明没有改进余地了,停止。
这种方法被称为**策略迭代**。
【算法 11.1】策略迭代
----
输入:策略 $\pi$,折扣 $\gamma$, 幕数 Episodes
1. 初始化
  初始化数组:$Q(S,A) \leftarrow 0$
  初始化策略:$\pi(S) \leftarrow Random$ 随机策略
2. 策略评估(多幕 Episodes 循环):
  根据当前策略与环境进行多幕交互
  得到 $Q(S,A)$
3. 策略改进
  $old.policy \leftarrow \pi$
  对每个状态更新策略:$\pi(s) \leftarrow \argmax_a Q(s)$
  如果 $\pi \ne old.policy$,跳转到第 2 步
  否则结束
输出:$Q(S,A), \pi$
----
在上述算法中,策略改进是我们本章要研究的重点,其它部分利用以前的知识都可以抽象成通用模块,所以我们来做一下模块划分,以便可以统一接口,专心研究算法部分。
模块可以做如图 11.1.1 的划分。
<center>
<img src="./img/PolicyIterationFlow.png">
图 11.1.1 策略迭代示意图
(左图:初始化在迭代中;右图:初始化在迭代外)
</center>
有左右两个子图,我们稍后再解释,先看看各个模块的具体功能。
#### 初始化模块
```python
# 策略迭代
class Policy_Iteration(object):
# 初始化
def __init__(self, env, policy, gamma, episodes, final, check_policy_method=0):
self.policy = policy # 初始策略
self.rough_episodes = episodes # 粗预测的分幕循环次数
self.final_episodes = final # 细预测的分幕循环次数
self.env = env # 环境
self.gamma = gamma # 折扣
self.check_method = check_policy_method # 迭代终止条件
self.nA = self.env.action_space.n # 动作空间
self.nS = self.env.observation_space.n # 状态空间
self.n_iteration = 0 # 分幕循环次数总和
self.Value = np.zeros((self.nS, self.nA)) # G 的总和
self.Count = np.zeros((self.nS, self.nA)) # G 的数量
```
初始化工作很简单,主要是接收初始化参数,其中最主要的工作是最后两行代码,把计算 G 值的数据清零。
#### 策略 $\pi$ 评估(粗预测)模块
本模块实际上就是每次访问法。代码与 10.5 节中的算法没有区别,就不再赘述了。
为什么叫做粗预测?在给定一个初始策略,或者是非最优策略后,我们要估算出它的 Q 值。由于蒙特卡洛方法需要很多次分幕采样才能得到较为精确的结果,所以耗时很长。但是,我们从 10.5 节中对于三个策略的预测可以看到,即使是很差的策略,在循环到一定的幕数后,也可以得到不错的方向性指导。这就提示我们,可以指定一个不是很大的循环次数,粗略得到 Q 表格即可。即使其中有些动作方向不是最优的,但是还有机会在后续的迭代中改进。
可以用一个简单的方格世界的例子来说明这个概念。
<center>
<img src="./img/GridWorld44.png">
图 11.1.2 简单的方格世界
</center>
在图 11.1.2 所示的 4x4 方格世界中,可以从任意点出发,终点为 $s_{15}$。在到达 $s_{15}$ 之前,每一步移动都得到 -1 的奖励,到 $s_{15}$ 的奖励为 0。智能体可以选择四个方向移动出界后会返回原地。无折扣。
由于存在 $s_{15}$ 终止状态,所以是一个分幕任务,可以用蒙特卡洛法来解决。当然,这个问题也可以用动态规划来解决,因为知道环境信息。
我们先看看用动态规划法得到的解:
<center>
<img src="./img/GridWorld44_DP.png">
图 11.1.2 简单的方格世界
Q 表格函数值;右:动作方向)
</center>
可以看到这种方法的结果是沿对角线成镜像对称的,很漂亮。
我们再看看用蒙特卡洛法循环 100幕、200幕、300幕、400幕、500幕和1000幕得到的解
<center>
<img src="./img/GridWorld44_MC.png">
图 11.1.2 简单的方格世界
</center>
- 100幕的结果只有左上角的状态位不靠谱
- 200幕时仍然是左上角的动作有误差
- 300幕时在细节上与后面的图还有差距但是在每个状态上的动作选择都是合理的尽管不是最优的
- 后续三图,动作选择已经没有什么可挑剔的了。
在比较图 11.1.2 和 图 11.1.3 后,读者可能会有疑问:在图 11.1.2 中,$s_0,s_5,s_{10}$ 三个状态上的动作是可以向下和向右,两个 q 函数值相等,而整体上看,也是以对角线为中心成镜像的。那为什么在图 11.1.3 中,即使是经过了 1000 幕,也没有出现这种完全对称的情况?
因为图 11.1.2 是用动态规划算法得到的结果,比较精准,如果一个状态中有两个动作等价,那这种方法确实可以得出两个相等的 q 值来。而图 11.1.3 是用蒙特卡洛法,就不那么精准了,误差在小数点以后两位都是正常的,所以极少能够得到两个完全相等的 q 值。
比如 $q_1$ = 0.0506$q_2$ = 0.0507,在我们看来 $q_1,q_2$ 之间的误差是蒙特卡洛方法造成的,可以忽略,但是算法就是认为它们不等,除非强制截断到小数点后两位。但是,不同的问题,或者相同的问题但折扣不同,小数点后截断多少位不是固定的,需要具体问题具体分析。
回到我们要解决的问题上来,事实证明,经过 300 幕的循环,就能够得到正确的策略了,在此之上继续做策略迭代没有问题,不必非得要循环到 Q 函数值收敛为止。而在实际的应用中,我们甚至可以让粗预测模块循环更少的次数,以加快整体迭代速度。
在图 11.1.1 所示的价值迭代的算法中(贝尔曼最优方程),每计算一次 Q 值,立刻从中取最大者来计算一次 V 值,相当于评估与改进的比例为 1:1这是一种极端的例子。
#### 策略改进模块
我们可以借助图 11.1.2 来理解策略改进机制。
<center>
<img src="./img/PolicyIteration.png">
图 11.1.2 策略迭代示意图
</center>
初始化策略可以命名为 $\pi_0$,而 $Q_0$ 在代码中虽然不存在,但是 $Q=Value/Count$,所以 self.Value=np.zeros(...) 和 self.Count=np.zeros(...) 代表了 $Q_0 \leftarrow 0$。
1. 初始化 $Q, \pi$ 为 0不妨命名为 $Q_0, \pi_0$
2. 使用蒙特卡洛每次访问法根据 $\pi_0$ 估算 $Q_1$
3. 使用贪心算法从 $Q_1$ 中抽取出新的策略 $\pi_1$
4. 比较 $\pi_0$ 和 $\pi_1$,如果一致则认为收敛;
5. 使用蒙特卡洛每次访问法根据 $\pi_1$ 估算 $Q_2$
6. 使用贪心算法从 $Q_2$ 中抽取出新的策略 $\pi_2$
7. 比较 $\pi_1$ 和 $\pi_2$,如果一致则认为收敛;
......
以此类推,直到收敛为止,就认为得到了最优策略 $\pi_*$,如果需要的话可以再估算一次 $Q$ 值来得到 $Q_*$
所以,所谓的策略改进模块(算法),就是 $\pi=greedy(Q)$ 贪心算法,也是我们本章研究的对象。
这里的代码只是定义了一个 policy_improvement() 函数,并没有具体实现,要靠后面将要研究的各种算法来重写此函数。
```python
# 策略改进(算法重写此函数)
def policy_improvement(self, Q):
pass
```
该函数输入为 $Q$ 表格值,实现时需要在内部修改 self.policy 的值,来改进策略。
#### 策略 $\pi_*$ 评估(细预测)模块
这一部分是可选的,用来验证我们认为的最优策略 $\pi_*$ 是不是名副其实。简单的验证方法是从 $\pi_*$ 得到 $Q_*$ 后,看看 $Q_*$ 的值是不是大于以往每个策略上得到的相应位置的 $Q$ 值,或者看循环多次以后,$Q_*$ 表格中的每个值是不是都可以收敛到某个具体的数值。
####
左图
初始化放在迭代内部,理由是:使用 $\pi$ 策略评估模块进行预测,生成 Q 表格后,通过策略改进模块抽取出新的策略 $\pi'$。此时应该再次初始化,以保证
右图
https://zhuanlan.zhihu.com/p/28498261
理论证明策略迭代方法的正确
- 可以比较policy iteration and value iteration
- 可以绘图:结果图、曲线(趋势)图
- maze https://cs.stanford.edu/people/karpathy/reinforcejs/gridworld_dp.html
$$
\pi_*=\argmax_\pi v_\pi(s)
$$

Просмотреть файл

@ -0,0 +1,26 @@
## 贪心法与探索性出发
#### 贪心法
#### 探索性出发Exploring Starts
在 8.4 节中,我们学习过 Q 函数的定义。如果应用到冰面行走问题上Q 函数表格会如表 10.4.1 所示。
表 10.4.1 二十一点问题的 Q 函数表格
|状态 $\downarrow$ 动作$\to$|HIT|Stick|
|:-:|:-:|:-:|
|12|$q_\pi(12,a_0)$|$q_\pi(12,a_1)$|$q_\pi(s_0,a_2)$|$q_\pi(s_0,a_3)$|
|13|$q_\pi(13,a_0)$|$q_\pi(13,a_1)$|$q_\pi(s_1,a_2)$|$q_\pi(s_1,a_3)$|
|...|...|...|...|...|
|20|$q_\pi(20,a_0)$|$q_\pi(20,a_1)$|
|21|$q_\pi(21,a_0)$|$q_\pi(21,a_1)$|
状态-动作的组合构成 Q 表格,一共有 10x2=20 个组合。如果想评价在某个状态上哪个动作最好,那么最起码要在该状态上尝试完所有动作后才会有评价的基础。

Просмотреть файл

@ -0,0 +1,11 @@
### 10.4.2 软性策略
$\varepsilon$-Soft
我们已经不是第一次接触软性策略了,在第二章的躲避强盗问题中,曾经学习过梯度上升法,里面使用了 Softmax 函数,根据各个动作的价值计算出备选概率,而不是使用非黑即白的硬性策略 argmax() 来选择后续动作。这样做的好处是:一方面在以最大概率选择(利用)了历史表现最好的动作的同时,给其它表现不好的动作一定的机会来进行探索。
### 10.4.3 GLIE 方法

Просмотреть файл

@ -0,0 +1,4 @@
on-policy
off-policy

Просмотреть файл

@ -0,0 +1,37 @@
- 每次访问,首次访问,估算 MRP 的 V
- MDP 平稳
- 试探性出发,可以被 e-soft 代替
- GLIM, e=1/k, k 为迭代次数k 越大信心越强e越小
- 非平稳
- 同轨策略控制
- 离轨策略 与 重要度采样
### 非平稳环境
什么是非平稳环境?
我们前面学习的例子,基本上都是平稳环境,即:到达状态 $s_t$ 后,一定可以得到奖励 $R_{t+1}$。下面我们举一些非平稳环境的例子。
- 良好路况下驾驶汽车,因为路面的摩擦力(可以认为是)恒定,司机可以按照通常的经验进行加速、转弯、刹车等动作。但是一旦遇到雨雪天气,以上动作都会导致车俩不会按照预期行驶,加速时前轮会空转,转弯时车辆会侧滑,刹车时需要更长的制动距离。
- 打游戏与程序Program对战时程序会按照预先设计好的一些套路来进攻或者防守摸清这些套路人类玩家就会胜利。但是如果与人工智能AI对战时对手会根据人类玩家的应对方式实时地修改策略这样前面学习的经验就没用了必须随机应变。
- 在股票市场中,一般的散户或者小机构可以按照自己的策略进行正常投资,然后根据市场行情进行调整、验证。但是大的机构会遇到一些麻烦,比如一个基金公司有 100 亿的规模,那么它的买卖行为会对某只股票、某个板块、甚至对大盘造成影响,它的策略必须考虑到这种影响。
以上这些例子,如果用到强化学习的环境中,都会造成当遇到 $s_t$ 状态时,有时会有很高的正向奖励,有时却是很低的正向奖励,甚至是负向奖励。
### 增量更新
### 批量更新算法
$$
V(s) = V(s) + \alpha (G_t - V)
$$

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 563 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 16 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 23 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 10 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 23 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 31 KiB

Двоичный файл не отображается.

После

Ширина:  |  Высота:  |  Размер: 18 KiB

Просмотреть файл

@ -0,0 +1,196 @@
import numpy as np
import tqdm
import math
# MC 策略评估(预测):每次访问法估算 V_pi
def MC_EveryVisit_V_Policy(env, episodes, gamma, policy):
nS = env.observation_space.n
nA = env.action_space.n
Value = np.zeros(nS) # G 的总和
Count = np.zeros(nS) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
Episode = [] # 一幕内的(状态,奖励)序列
s = env.reset() # 重置环境,开始新的一幕采样
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, _ = env.step(action)
Episode.append((s, reward))
s = next_s
# 从后向前遍历计算 G 值
G = 0
for t in range(len(Episode)-1, -1, -1):
s, r = Episode[t]
G = gamma * G + r
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
# 检查是否收敛
Count[Count==0] = 1 # 把分母为0的填成1主要是对终止状态
V = Value / Count # 求均值
return V
# MC 策略评估(预测):每次访问法估算 Q_pi
def MC_EveryVisit_Q_Policy(env, episodes, gamma, policy):
nA = env.action_space.n # 动作空间
nS = env.observation_space.n # 状态空间
Value = np.zeros((nS, nA)) # G 的总和
Count = np.zeros((nS, nA)) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
# 重置环境,开始新的一幕采样
s = env.reset()
Episode = [] # 一幕内的(状态,动作,奖励)序列
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, _ = env.step(action)
Episode.append((s, action, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s,a] += G # 值累加
Count[s,a] += 1 # 数量加 1
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
Q = Value / Count # 求均值
return Q
def get_start_state(env, start_state):
s = env.reset()
while s != start_state: # 以随机选择的 s 为起始状态
action = np.random.choice(env.action_space.n)
next_s, _, done, _ = env.step(action)
s = next_s
if s == start_state:
break
if done is True:
s = env.reset()
return s
# MC 控制 探索出发
def MC_ES(env, episodes, gamma, policy):
nA = env.action_space.n
nS = env.observation_space.n
Value = np.zeros((nS, nA)) # G 的总和
Count = np.zeros((nS, nA)) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
# 重置环境,开始新的一幕采样
start_state = np.random.choice(nS)
s = get_start_state(env, start_state)
assert(s == start_state)
# 找到了指定的 start_state开始采样
Trajectory = [] # 一幕内的(状态,奖励)序列
action = np.random.choice(nA) # 起始动作也是随机的
next_s, reward, done, _ = env.step(action)
Trajectory.append((s, action, reward))
s = next_s
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s]) # 根据改进的策略采样
next_s, reward, done, _ = env.step(action)
Trajectory.append((s, action, reward))
s = next_s
num_step = len(Trajectory)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Trajectory[t]
G = gamma * G + r
Value[s,a] += G # 值累加
Count[s,a] += 1 # 数量加 1
Count[Count==0] = 1
Q = Value / Count # 求均值
for s in range(nS):
max_A = np.max(Q[s])
argmax_A = np.where(Q[s] == max_A)[0]
A = np.random.choice(argmax_A)
policy[s] = 0
policy[s,A] = 1
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
Q = Value / Count # 求均值
return Q
# GLIE - Greedy in the Limit with Infinite Exploration
def MC_EveryVisit_Q_GLIE(env, start_state, episodes, gamma, policy):
nA = env.action_space.n
nS = env.observation_space.n
Q = np.zeros((nS, nA)) # Q 动作价值
N = np.zeros((nS, nA)) # N 次数
for k in tqdm.trange(episodes): # 多幕循环
# 重置环境,开始新的一幕采样
s, info = env.reset(return_info=True)
Trajectory = [] # 一幕内的(状态,奖励)序列
s = start_state
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, info = env.step(action)
Trajectory.append((s, action, reward))
s = next_s
num_step = len(Trajectory)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Trajectory[t]
G = gamma * G + r
N[s,a] += 1 # 数量加 1
Q[s,a] += (G - Q[s,a])/N[s,a]
epsilon = 1 / (math.log(k+1)+1)
other_p = epsilon / nA
best_p = 1 - epsilon + epsilon/nA
# 更新策略
for s in range(nS):
max_A = np.max(Q[s])
argmax_A = np.where(Q[s] == max_A)[0]
A = np.random.choice(argmax_A)
policy[s] = other_p
policy[s,A] = best_p
return Q
# MC3 - 增量更新
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
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
# calculate G_t
G = 0
# 从后向前遍历
for t in range(len(trajectory)-1, -1, -1):
S, R = trajectory[t]
G = gamma * G + R
V[S] = V[S] + alpha * (G - V[S])
#endfor
#endfor
return V

Просмотреть файл

@ -0,0 +1,57 @@
import numpy as np
import copy
# 式 (9.6.1) 计算 q假设已知 v*(s')
def q_star(p_s_r_d, gamma, V):
q = 0
# 遍历每个转移概率,以计算 q
for p, s_next, reward, done in p_s_r_d:
# math: \sum_{s'} p_{ss'}^a [ r_{ss'}^a + \gamma * v_{\pi}(s')]
q += p * (reward + gamma * V[s_next] * (1-done))
return q
# 式 (9.6.3) 计算 v*
def v_star(s, actions, gamma, V, Q):
list_q = [] # 准备列表记录所有下游的 q*
for action in list(actions.keys()): # actions 是一个字典数据key 是动作
# 遍历每个动作以计算q值进而计算v值
q = q_star(actions[action], gamma, V) # 计算 q*
list_q.append(q) # 加入列表
Q[s,action] = q # 记录下所有的q(s,a)值,不需要再单独计算一次
return max(list_q) # 返回几个q*中的最大值,即 v=max(q)
def calculate_VQ_star(env, gamma, max_iteration):
V = np.zeros(env.observation_space.n) # 初始化 V(s)
Q = np.zeros((env.observation_space.n, env.action_space.n)) # 初始化 Q(s,a)
count = 0
# 迭代
while (count < max_iteration):
V_old = V.copy()
# 遍历所有状态 s
for s in range(env.observation_space.n):
actions = env.P[s]
V[s] = v_star(s, actions, gamma, V, Q)
# 检查收敛性
if abs(V-V_old).max() < 1e-4:
break
count += 1
# end while
print("迭代次数 = ",count)
return V, Q
def get_policy(env, V, gamma):
policy = np.zeros((env.nS, env.nA))
for s in range(env.nS):
actions = env.get_actions(s)
list_q = []
if actions is None:
continue
# 遍历每个策略概率
for action, next_p_s_r in actions:
q_star = 0
for p, s_next, r in next_p_s_r:
q_star += p * (r + gamma * V[s_next])
list_q.append(q_star)
policy[s, np.argmax(list_q)] = 1
return policy

Просмотреть файл

@ -0,0 +1,41 @@
import numpy as np
# 式 (8.4.2) 计算 q_pi
def q_pi(p_s_r_d, gamma, V):
q = 0
# 遍历每个转移概率,以计算 q_pi
for p, s_next, reward, done in p_s_r_d:
# math: \sum_{s'} p_{ss'}^a [ r_{ss'}^a + \gamma * v_{\pi}(s')]
q += p * (reward + gamma * V[s_next])
return q
# 式 (8.4.5) 计算 v_pi
def v_pi(policy, s, actions, gamma, V, Q):
v = 0
for action in list(actions.keys()): # actions 是一个字典数据key 是动作
q = q_pi(actions[action], gamma, V)
# math: \sum_a \pi(a|s) q_\pi (s,a)
v += policy[s][action] * q
# 顺便记录下q(s,a)值,不需要再单独计算一次
Q[s,action] = q
return v
# 迭代法计算 v_pi
def calculate_VQ_pi(env, policy, gamma, iteration, delta=1e-4):
V = np.zeros(env.observation_space.n) # 初始化 V(s)
Q = np.zeros((env.observation_space.n, env.action_space.n)) # 初始化 Q(s,a)
count = 0 # 计数器,用于衡量性能和避免无限循环
# 迭代
while (count < iteration):
V_old = V.copy() # 保存上一次的值以便检查收敛性
# 遍历所有状态 s
for s in range(env.observation_space.n):
actions = env.P[s]
V[s] = v_pi(policy, s, actions, gamma, V, Q)
# 检查收敛性
if abs(V-V_old).max() < delta:
break
count += 1
# end while
print("迭代次数 = ",count)
return V, Q

Просмотреть файл

@ -0,0 +1,82 @@
import numpy as np
import tqdm
# 策略迭代
class Policy_Iteration(object):
# 初始化
def __init__(self, env, policy_args, gamma, rough, final, check_policy_method=0):
self.rough_episodes = rough # 粗预测的分幕循环次数
self.final_episodes = final # 细预测的分幕循环次数
self.env = env # 环境
self.gamma = gamma # 折扣
self.check_method = check_policy_method # 迭代终止条件
self.nA = self.env.action_space.n # 动作空间
self.n_iteration = 0 # 分幕循环次数总和
if env.spec.id == 'Blackjack-v1':
# 0-21=22, 0-10=11, A/noA=2, Hit/Stick=2
self.Value = np.zeros((22, 11, 2, self.nA)) # G 的总和
self.Count = np.zeros((22, 11, 2, self.nA)) # G 的数量
else:
self.nS = self.env.observation_space.n # 状态空间
self.Value = np.zeros((self.nS, self.nA))
self.Count = np.zeros((self.nS, self.nA))
self.policy = self.init_policy(policy_args)
def init_policy(self, policy_args):
if self.env.spec.id == 'Blackjack-v1':
policy = np.zeros_like(self.Value)
policy[:,:,:] = policy_args
return policy
# 策略评估
def policy_evaluation(self, episodes):
for _ in tqdm.trange(episodes): # 多幕循环
self.n_iteration += 1
# 重置环境,开始新的一幕采样
s = self.env.reset()
Episode = [] # 一幕内的(状态,奖励)序列
done = False
while (done is False): # 幕内循环
int_s = (s[0], s[1], int(s[2]))
action = np.random.choice(self.nA, p=self.policy[int_s])
next_s, reward, done, _ = self.env.step(action)
Episode.append((int_s, action, reward))
s = next_s
# 从后向前遍历计算 G 值
G = 0
for t in range(len(Episode)-1, -1, -1):
s, a, r = Episode[t]
G = self.gamma * G + r
self.Value[s][a] += G # 值累加
self.Count[s][a] += 1 # 数量加 1
# 多幕循环结束,计算 Q 值
self.Count[self.Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
Q = self.Value / self.Count # 求均值
return Q
# 策略改进(算法重写此函数)
def policy_improvement(self, Q):
pass
# 策略迭代
def policy_iteration(self):
while True:
print("策略评估")
Q = self.policy_evaluation(self.rough_episodes)
print("策略改进")
old_policy = self.policy.copy()
self.policy_improvement(Q)
if self.check_policy_diff(old_policy, self.policy):
print("新旧策略相等")
break
print("精准的策略评估")
Q = self.policy_evaluation(self.final_episodes)
return Q, self.policy
def check_policy_diff(self, old_policy, new_policy):
if (self.check_method == 0):
return (old_policy == new_policy).all()
elif (self.check_method == 1):
old_arg = np.argmax(old_policy, axis=1)
new_arg = np.argmax(new_policy, axis=1)
return (old_arg == new_arg).all()

Просмотреть файл

@ -0,0 +1,104 @@
import numpy as np
import common.GridWorld_Model as model
import Algorithm.Algo_PolicyValueFunction as algoPi
import common.CommonHelper as helper
import common.DrawQpi as drawQ
import tqdm
# 状态空间 = 空间宽度 x 空间高度
GridWidth, GridHeight = 4, 4
# 起点,可以多个
StartStates = [0]
# 终点,可以多个
EndStates = [15]
# 动作空间
LEFT, DOWN, RIGHT, UP = 0, 1, 2, 3
Actions = [LEFT, DOWN, RIGHT, UP]
# 初始策略
Policy = [0.25, 0.25, 0.25, 0.25]
# 转移概率: [SlipLeft, MoveFront, SlipRight, SlipBack]
Transition = [0.0, 1.0, 0.0, 0.0]
# 每走一步的奖励值可以是0或者-1
StepReward = -1
# 特殊奖励 from s->s' then get r, 其中 s,s' 为状态序号,不是坐标位置
SpecialReward = {
#(1,0):0,
#(4,0):0,
(14,15):0,
(11,15):0
}
# 特殊移动,用于处理类似虫洞场景
SpecialMove = {
}
# 墙
Blocks = []
# MC 策略评估(预测):每次访问法估算 Q_pi
def MC_EveryVisit_Q_Policy_test(env, episodes, gamma, policy, checkpoint):
Q_history = []
nA = env.action_space.n # 动作空间
nS = env.observation_space.n # 状态空间
Value = np.zeros((nS, nA)) # G 的总和
Count = np.zeros((nS, nA)) # G 的数量
for episode in tqdm.trange(episodes): # 多幕循环
# 重置环境,开始新的一幕采样
s = env.reset()
Episode = [] # 一幕内的(状态,动作,奖励)序列
done = False
while (done is False): # 幕内循环
action = np.random.choice(nA, p=policy[s])
next_s, reward, done, _ = env.step(action)
Episode.append((s, action, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s,a] += G # 值累加
Count[s,a] += 1 # 数量加 1
if (episode + 1) in checkpoint:
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
Q = Value / Count # 求均值
Q = np.round(Q, 0)
Q_history.append(Q)
return Q_history
if __name__=="__main__":
np.random.seed(15)
env = model.GridWorld(
GridWidth, GridHeight, StartStates, EndStates, # 关于状态的参数
Actions, Policy, Transition, # 关于动作的参数
StepReward, SpecialReward, # 关于奖励的参数
SpecialMove, Blocks) # 关于移动的限制
#model.print_P(env.P_S_R)
policy = helper.create_policy(env.nS, env.nA, (0.25, 0.25, 0.25, 0.25))
gamma = 1
max_iteration = 1000
V_pi, Q_pi = algoPi.calculate_VQ_pi(env, policy, gamma, max_iteration, delta=1e-6)
helper.print_seperator_line(helper.SeperatorLines.long, "动态规划解")
helper.print_seperator_line(helper.SeperatorLines.short, "V 函数")
print(np.around(V_pi,0).reshape(4,4))
helper.print_seperator_line(helper.SeperatorLines.short, "Q 函数")
print(np.around(Q_pi,0))
drawQ.draw(Q_pi, (4,4))
checkpoint = [100,200,300,400,500,1000]
policy = helper.create_policy(env.nS, env.nA, (0.25, 0.25, 0.25, 0.25))
Q_history = MC_EveryVisit_Q_Policy_test(env, max_iteration, gamma, policy, checkpoint)
assert(len(checkpoint)==len(Q_history))
for i in range(len(checkpoint)):
helper.print_seperator_line(helper.SeperatorLines.long, str(checkpoint[i]))
Q = Q_history[i]
V = helper.calculat_V_from_Q(Q, policy)
helper.print_seperator_line(helper.SeperatorLines.short, "V 函数")
print(np.round(V,0).reshape(4,4))
helper.print_seperator_line(helper.SeperatorLines.short, "Q 函数")
print(np.around(Q, 0))
drawQ.draw(Q, (4,4))

Просмотреть файл

@ -0,0 +1,93 @@
import numpy as np
import common.GridWorld_Model as model
import common.CommonHelper as helper
import common.DrawQpi as drawQ
import Algorithm.Base_MC_Policy_Iteration as algoMC
import Algorithm.Algo_MC_Policy_Evaulation as algoV
# 状态空间 = 空间宽度 x 空间高度
GridWidth, GridHeight = 4, 4
# 起点,可以多个
StartStates = [5,6,9,10]
# 终点,可以多个
EndStates = [0,15]
# 动作空间
LEFT, DOWN, RIGHT, UP = 0, 1, 2, 3
Actions = [LEFT, DOWN, RIGHT, UP]
# 初始策略
Policy = [0.25, 0.25, 0.25, 0.25]
# 转移概率: [SlipLeft, MoveFront, SlipRight, SlipBack]
Transition = [0.0, 1.0, 0.0, 0.0]
# 每走一步的奖励值可以是0或者-1
StepReward = -1
# 特殊奖励 from s->s' then get r, 其中 s,s' 为状态序号,不是坐标位置
SpecialReward = {
(1,0):0,
(4,0):0,
(14,15):0,
(11,15):0
}
# 特殊移动,用于处理类似虫洞场景
SpecialMove = {
}
# 墙
Blocks = []
class MC_E_Greedy(algoMC.Policy_Iteration):
def __init__(self, env, policy, gamma:float, rough:int, final:int, epsilon:float):
super().__init__(env, policy, gamma, rough, final)
self.epsilon = epsilon
self.other_p = self.epsilon / self.nA
self.best_p = 1 - self.epsilon + self.epsilon / self.nA
def policy_improvement(self, Q):
print(np.sum(Q))
for s in range(self.nS):
if s in env.EndStates:
self.policy[s] = 0
else:
max_A = np.max(Q[s])
argmax_A = np.where(Q[s] == max_A)[0]
A = np.random.choice(argmax_A)
self.policy[s] = self.other_p
self.policy[s,A] = self.best_p
return self.policy
if __name__=="__main__":
gamma = 1
rough_episodes = 500
final_episodes = 10000
epsilon = 0.1
np.random.seed(15)
env = model.GridWorld(
GridWidth, GridHeight, StartStates, EndStates, # 关于状态的参数
Actions, Policy, Transition, # 关于动作的参数
StepReward, SpecialReward, # 关于奖励的参数
SpecialMove, Blocks) # 关于移动的限制
policy = helper.create_policy(env.nS, env.nA, (0.25,0.25,0.25,0.25))
policy_0 = policy.copy()
env.reset()
algo = MC_E_Greedy(env, policy, gamma, rough_episodes, final_episodes, epsilon)
Q, policy = algo.policy_iteration()
env.close()
print("------ 最优动作价值函数 -----")
Q=np.round(Q,0)
print(Q)
drawQ.draw(Q,(4,4))
print(policy)
env.reset()
V = algoV.MC_EveryVisit_V_Policy(env, 10000, gamma, policy_0)
print(np.reshape(V, (4,4)))
env.reset()
V = algoV.MC_EveryVisit_V_Policy(env, 10000, gamma, policy)
print(np.reshape(V, (4,4)))
V = helper.calculat_V_from_Q(Q, policy)
print(np.reshape(V, (4,4)))

Просмотреть файл

@ -0,0 +1,53 @@
import numpy as np
import gym
import Algorithm.Algo_MC_Policy_Iteration as algoMC
import Algorithm.Algo_OptimalValueFunction as algoDP
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
class MC_Greedy(algoMC.Policy_Iteration):
def policy_improvement(self, Q):
for s in range(self.nS):
arg = np.argmax(Q[s])
self.policy[s] = 0
self.policy[s, arg] = 1
return self.policy
def get_groud_truth(env, gamma):
iteration = 100
_, Q = algoDP.calculate_VQ_star(env, gamma, iteration)
return Q
if __name__=="__main__":
gamma = 0.9
episodes = 1000
final = 2000
np.set_printoptions(suppress=True)
env = gym.make("FrozenLake-v1", desc=None, map_name = "8x8", is_slippery=False)
Q_real = get_groud_truth(env, gamma)
print(np.round(Q_real,3))
drawQ.draw(Q_real,(8,8))
policy = helper.create_policy(env, (0.25,0.25,0.25,0.25))
env.reset(seed=5)
algo = MC_Greedy(env, policy, gamma, episodes, final)
Q, policy = algo.policy_iteration()
env.close()
print("------ 最优动作价值函数 -----")
print(np.round(Q,3))
drawQ.draw(policy,(8,8))
print(policy)

Просмотреть файл

@ -0,0 +1,70 @@
from cProfile import label
from email import policy
import numpy as np
import gym
import Algorithm.Base_MC_Policy_Iteration as algoMC
import Algorithm.Algo_OptimalValueFunction as algoDP
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
class MC_E_Greedy(algoMC.Policy_Iteration):
def __init__(self, env, policy, gamma:float, episodes:int, final:int, epsilon:float):
super().__init__(env, policy, gamma, episodes, final)
self.epsilon = epsilon
self.other_p = self.epsilon / self.nA
self.best_p = 1 - self.epsilon + self.epsilon / self.nA
def policy_improvement(self, Q):
print(np.sum(Q))
for s in range(self.nS):
if s in end_states:
self.policy[s] = 0
else:
max_A = np.max(Q[s])
argmax_A = np.where(Q[s] == max_A)[0]
A = np.random.choice(argmax_A)
self.policy[s] = self.other_p
self.policy[s,A] = self.best_p
return self.policy
def get_groud_truth(env, gamma):
iteration = 100
_, Q = algoDP.calculate_VQ_star(env, gamma, iteration)
return Q
if __name__=="__main__":
gamma = 0.9
episodes = 10
final = 2000
epsilon = 0.05
np.set_printoptions(suppress=True)
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=False)
end_states = [5, 7, 11, 12, 15]
Q_real = get_groud_truth(env, gamma)
print(np.round(Q_real,3))
drawQ.draw(Q_real,(4,4))
policy = helper.create_policy(env.observation_space.n, env.action_space.n, (0.25,0.25,0.25,0.25))
env.reset(seed=5)
algo = MC_E_Greedy(env, policy, gamma, episodes, final, epsilon)
Q, policy = algo.policy_iteration()
env.close()
print("------ 最优动作价值函数 -----")
print(np.round(Q,3))
drawQ.draw(Q,(4,4))
print(policy)

Просмотреть файл

@ -0,0 +1,67 @@
from cProfile import label
from email import policy
import numpy as np
import gym
import Algorithm.Algo_MC_Policy_Iteration as algoMC
import Algorithm.Algo_OptimalValueFunction as algoDP
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
class MC_GLIE(algoMC.Policy_Iteration):
def policy_improvement(self, Q):
print(np.sum(Q))
epsilon = 1 / (math.log(self.n_iteration+1)+1)
other_p = epsilon / self.nA
best_p = 1 - epsilon + epsilon/self.nA
for s in range(self.nS):
if s in end_states:
self.policy[s] = 0
else:
max_A = np.max(Q[s])
argmax_A = np.where(Q[s] == max_A)[0]
A = np.random.choice(argmax_A)
self.policy[s] = other_p
self.policy[s,A] = best_p
return self.policy
def get_groud_truth(env, gamma):
iteration = 100
_, Q = algoDP.calculate_VQ_star(env, gamma, iteration)
return Q
if __name__=="__main__":
gamma = 0.9
episodes = 10
final = 2000
np.set_printoptions(suppress=True)
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=False)
end_states = [5, 7, 11, 12, 15]
Q_real = get_groud_truth(env, gamma)
print(np.round(Q_real,3))
drawQ.draw(Q_real,(4,4))
policy = helper.create_policy(env, (0.25,0.25,0.25,0.25))
env.reset(seed=5)
algo = MC_GLIE(env, policy, gamma, episodes, final, check_policy_method=1)
Q, policy = algo.policy_iteration()
env.close()
print("------ 最优动作价值函数 -----")
print(np.round(Q,3))
drawQ.draw(Q,(4,4))
print(policy)

Просмотреть файл

@ -0,0 +1,109 @@
from cProfile import label
from email import policy
import numpy as np
import gym
import Algorithm.Algo_MC_Policy_Iteration as algoMC
import Algorithm.Algo_OptimalValueFunction as algoDP
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import tqdm
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
class MC_Greedy(algoMC.Policy_Iteration):
def policy_improvement(self, Q):
for s in range(self.nS):
arg = np.argmax(Q[s])
self.policy[s] = 0
self.policy[s, arg] = 1
return self.policy
class MC_E_Greedy(algoMC.Policy_Iteration):
def __init__(self, env, policy, episodes, gamma, epsilon):
super().__init__(env, policy, episodes, gamma)
self.epsilon = epsilon
def initialize(self):
super().initialize()
self.other_p = self.epsilon / self.nA
self.best_p = 1 - self.epsilon + self.epsilon / self.nA
def policy_improvement(self, Q):
for s in range(self.nS):
max_A = np.max(Q[s])
if max_A == 0:
self.policy[s] = 0
else:
argmax_A = np.where(Q[s] == max_A)[0]
A = np.random.choice(argmax_A)
self.policy[s] = self.other_p
self.policy[s,A] = self.best_p
return self.policy
class MC_GLIE(algoMC.Policy_Iteration):
def __init__(self, env, policy, episodes, gamma, epsilon):
super().__init__(env, policy, episodes, gamma)
self.epsilon = epsilon
def initialize(self):
super().initialize()
epsilon = 1
self.other_p = epsilon / self.nA
self.best_p = 1 - epsilon + epsilon/self.nA
def policy_improvement(self, Q):
epsilon = 1 / (math.log(k+1)+1)
other_p = epsilon / self.nA
best_p = 1 - epsilon + epsilon/self.nA
# 更新策略
for s in range(self.nS):
max_A = np.max(Q[s])
if max_A == 0:
self.policy[s] = 0
else:
argmax_A = np.where(Q[s] == max_A)[0]
A = np.random.choice(argmax_A)
self.policy[s] = self.other_p
self.policy[s,A] = self.best_p
return self.policy
def get_groud_truth(env, gamma):
iteration = 100
_, Q = algoDP.calculate_VQ_star(env, gamma, iteration)
return Q
if __name__=="__main__":
gamma = 0.9
episodes = 1000
np.set_printoptions(suppress=True)
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=False)
Q_real = get_groud_truth(env, gamma)
print(np.round(Q_real,3))
drawQ.draw(Q_real,(4,4))
policy = helper.create_policy(env, (0.25,0.25,0.25,0.25))
env.reset(seed=5)
algo = MC_Greedy(env, policy, episodes, gamma)
#algo = MC_E_Greedy(env, policy, episodes, gamma, 0.1)
Q, policy = algo.policy_iteration()
env.close()
print("------ 最优动作价值函数 -----")
print(np.round(Q,3))
drawQ.draw(Q,(4,4))
print(policy)

Просмотреть файл

@ -0,0 +1,19 @@
import gym
import numpy as np
env = gym.make('Blackjack-v1', sab=True)
print(env.action_space)
print(env.observation_space)
for i in range(100):
s = env.reset()
Episode = []
while True:
action = np.random.choice(2)
next_s, r, done, _ = env.step(action)
Episode.append((s, action, r, next_s))
if done:
break
s = next_s
print(Episode)
env.close()

Просмотреть файл

@ -0,0 +1,71 @@
import numpy as np
import gym
import Algorithm.Algo_OptimalValueFunction as algoDP
import Algorithm.Base_MC_Policy_Iteration as base
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
class MC_E_Greedy(base.Policy_Iteration):
def __init__(self, env, policy, gamma:float, episodes:int, final:int, epsilon:float):
super().__init__(env, policy, gamma, episodes, final)
self.epsilon = epsilon
self.other_p = self.epsilon / self.nA
self.best_p = 1 - self.epsilon + self.epsilon / self.nA
def policy_improvement(self, Q):
print(np.sum(Q))
for s in range(self.nS):
if s in end_states:
self.policy[s] = 0
else:
max_A = np.max(Q[s])
argmax_A = np.where(Q[s] == max_A)[0]
A = np.random.choice(argmax_A)
self.policy[s] = self.other_p
self.policy[s,A] = self.best_p
return self.policy
def get_groud_truth(env, gamma):
iteration = 100
_, Q = algoDP.calculate_VQ_star(env, gamma, iteration)
return Q
if __name__=="__main__":
gamma = 0.9
rough = 300
final = 10000
epsilon = 0.05
np.set_printoptions(suppress=True)
#desc = ["SFHF","FFFF","HFHF","FFFF","FFFG"]
env = gym.make("FrozenLake-v1", desc=None, map_name = "4x4", is_slippery=True)
end_states = [5, 7, 11, 12, 15]
helper.print_seperator_line(helper.SeperatorLines.long, "动态规划法")
Q_real = get_groud_truth(env, gamma)
helper.print_seperator_line(helper.SeperatorLines.short, "Q 函数")
print(np.round(Q_real,3))
drawQ.draw(Q_real,(4,4))
policy = helper.create_policy(env.observation_space.n, env.action_space.n, (0.25,0.25,0.25,0.25))
env.reset(seed=5)
algo = MC_E_Greedy(env, policy, gamma, rough, final, epsilon)
Q, policy = algo.policy_iteration()
env.close()
helper.print_seperator_line(helper.SeperatorLines.long, "蒙特卡洛法 - E-Greedy")
helper.print_seperator_line(helper.SeperatorLines.short, "Q 函数")
print(np.round(Q,3))
drawQ.draw(Q,(4,4))
print(policy)

Просмотреть файл

@ -0,0 +1,45 @@
from email import policy
import numpy as np
import gym
import Algorithm.Algo_OptimalValueFunction as algoDP
import Algorithm.Base_MC_Policy_Iteration as base
import common.DrawQpi as drawQ
import common.CommonHelper as helper
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
class MC_Greedy(base.Policy_Iteration):
def policy_improvement(self, Q):
for i in range(Q.shape[0]):
for j in range(Q.shape[1]):
for k in range(Q.shape[2]):
arg = np.argmax(Q[i,j,k])
self.policy[i,j,k] = 0
self.policy[i,j,k,arg] = 1
return self.policy
def get_groud_truth(env, gamma):
iteration = 100
_, Q = algoDP.calculate_VQ_star(env, gamma, iteration)
return Q
if __name__=="__main__":
gamma = 1
rough = 1000
final = 20000
np.set_printoptions(suppress=True)
env = gym.make('Blackjack-v1', sab=True)
env.reset(seed=5)
algo = MC_Greedy(env, (0.5,0.5), gamma, rough, final)
Q, policy = algo.policy_iteration()
env.close()
print(policy)

Просмотреть файл

@ -0,0 +1,137 @@
import numpy as np
import gym
import Algorithm.Base_MC_Policy_Iteration as base
import matplotlib.pyplot as plt
import seaborn as sns
import tqdm
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
HIT = 1
STICK = 0
def static_policy(state):
if state[0] >= 20:
return STICK
else:
return HIT
def run_V(env, gamma, episodes):
# 玩家手牌点数和1-21, 单元 0 空置
# 庄家明牌点数: 1-10, 单元 0 空置
# 有无可用的A 0-1 (无/有)
Value = np.zeros((22, 11, 2)) # G 的总和
Count = np.zeros((22, 11, 2)) # G 的数量
for episode in tqdm.trange(episodes):
s = env.reset()
Episode = [] # 一幕内的(状态,动作,奖励)序列
done = False
while (done is False): # 幕内循环
int_s = (s[0], s[1], int(s[2]))
action = static_policy(s)
next_s, reward, done, _ = env.step(action)
Episode.append((int_s, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, r = Episode[t]
G = gamma * G + r
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
print(Count)
print(np.min(Count[12:,1:,]))
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
V = Value / Count # 求均值
return V
def greedy_policy(state, Q):
max_A = np.max(Q[state]) # 从几个动作中取最大值
argmax_A = np.where(Q[state] == max_A)[0]
action = np.random.choice(argmax_A)
return action
def run_Q(env, gamma, episodes, policy):
# 玩家手牌点数和1-21, 单元 0 空置
# 庄家明牌点数: 1-10, 单元 0 空置
# 有无可用的A 0-1 (无/有)
# 动作0-停牌1-要牌
Value = np.zeros((22, 11, 2, 2)) # G 的总和
Count = np.ones((22, 11, 2, 2)) # G 的数量
Q = Value / Count
for episode in tqdm.trange(episodes):
s = env.reset()
Episode = [] # 一幕内的(状态,动作,奖励)序列
done = False
while (done is False): # 幕内循环
int_s = (s[0], s[1], int(s[2]))
action = greedy_policy(s, Q)
next_s, reward, done, _ = env.step(action)
Episode.append((int_s, action, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s][a] += G # 值累加
Count[s][a] += 1 # 数量加 1
# refine policy
max_A = np.max(Q[state]) # 从几个动作中取最大值
argmax_A = np.where(Q[state] == max_A)[0]
action = np.random.choice(argmax_A)
print(Count)
print(np.min(Count[12:,1:,]))
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
V = Value / Count # 求均值
return V
if __name__=="__main__":
gamma = 1
episodes = 50000
env = gym.make('Blackjack-v1', sab=True)
V = run_V(env, gamma, episodes)
env.close()
#print(V)
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(121)
sns.heatmap(np.flipud(V[12:,1:,0]), cmap="Wistia", xticklabels=range(1, 11),
yticklabels=list(reversed(range(12, 22))))
ax.set_ylabel('玩家的手牌', fontsize=16)
ax.set_xlabel('庄家的明牌', fontsize=16)
ax.set_title('无可用的 A')
ax = fig.add_subplot(122)
sns.heatmap(np.flipud(V[12:,1:,1]), cmap="Wistia", xticklabels=range(1, 11),
yticklabels=list(reversed(range(12, 22))))
ax.set_ylabel('玩家的手牌', fontsize=16)
ax.set_xlabel('庄家的明牌', fontsize=16)
ax.set_title('有可用的 A')
plt.show()
fig = plt.figure()
ax = Axes3D(fig)
ax.set_zticks(np.arange(-1,1,1))
X, Y = np.meshgrid(range(0,10), range(0,10))
ax.plot_surface(X, Y, V[12:,1:,0])
plt.show()

Просмотреть файл

@ -0,0 +1,182 @@
import numpy as np
import gym
import Algorithm.Base_MC_Policy_Iteration as base
import matplotlib.pyplot as plt
import seaborn as sns
import tqdm
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
HIT = 1
STICK = 0
def static_policy(state):
if state[0] >= 20:
return STICK
else:
return HIT
def run_V(env, gamma, episodes):
# 玩家手牌点数和1-21, 单元 0 空置
# 庄家明牌点数: 1-10, 单元 0 空置
# 有无可用的A 0-1 (无/有)
Value = np.zeros((22, 11, 2)) # G 的总和
Count = np.zeros((22, 11, 2)) # G 的数量
for episode in tqdm.trange(episodes):
s = env.reset()
Episode = [] # 一幕内的(状态,动作,奖励)序列
done = False
while (done is False): # 幕内循环
int_s = (s[0], s[1], int(s[2]))
action = static_policy(s)
next_s, reward, done, _ = env.step(action)
Episode.append((int_s, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, r = Episode[t]
G = gamma * G + r
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
print(Count)
print(np.min(Count[12:,1:,]))
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
V = Value / Count # 求均值
return V
def greedy_policy(state, Q):
max_A = np.max(Q[state]) # 从几个动作中取最大值
argmax_A = np.where(Q[state] == max_A)[0]
action = np.random.choice(argmax_A)
return action
def run_Q_ES(env, gamma, episodes, policy):
# 玩家手牌点数和1-21, 单元 0 空置
# 庄家明牌点数: 1-10, 单元 0 空置
# 有无可用的A 0-1 (无/有)
# 动作0-停牌1-要牌
Value = np.zeros((22, 11, 2, 2)) # G 的总和
Count = np.ones((22, 11, 2, 2)) # G 的数量
for episode in tqdm.trange(episodes):
s = env.reset()
Episode = [] # 一幕内的(状态,动作,奖励)序列
# ES - initial random choice
int_s = (s[0], s[1], int(s[2]))
action = np.random.choice(2) # 这里必须随机取不能使用policy才符合探索的含义
next_s, reward, done, _ = env.step(action)
Episode.append((int_s, action, reward))
s = next_s # 迭代
while (done is False): # 幕内循环
int_s = (s[0], s[1], int(s[2])) # 把返回的 False/True 变成 0/1
action = np.random.choice(2, p=policy[int_s])
next_s, reward, done, _ = env.step(action)
Episode.append((int_s, action, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s][a] += G # 值累加
Count[s][a] += 1 # 数量加 1
# greedy policy
qs = Value[s] / Count[s]
policy[s] = 0
argmax = np.argmax(qs)
policy[s][argmax] = 1
return Value/Count
def create_blackjack_policy(args):
assert(args[0]+args[1]==1)
policy= np.zeros((22,11,2,2))
policy[:,:,:,0] = args[0]
policy[:,:,:,1] = args[1]
return policy
if __name__=="__main__":
policy = create_blackjack_policy((0.5,0.5))
gamma = 1
episodes = 500000
env = gym.make('Blackjack-v1', sab=True)
Q = run_Q_ES(env, gamma, episodes, policy)
env.close()
state_value_no_usable_ace = np.max(Q[12:, 1:, 0, :], axis=-1)
state_value_usable_ace = np.max(Q[12:, 1:, 1, :], axis=-1)
# get the optimal policy
action_no_usable_ace = np.argmax(Q[12:, 1:, 0, :], axis=-1)
action_usable_ace = np.argmax(Q[12:, 1:, 1, :], axis=-1)
print(Q[18,:,1,:])
images = [action_usable_ace,
state_value_usable_ace,
action_no_usable_ace,
state_value_no_usable_ace]
titles = ['最优策略(有可用的 A',
'最优状态价值(有可用的 A',
'最优策略(没有可用的 A',
'最优状态价值(没有可用的 A']
_, axes = plt.subplots(2, 2, figsize=(15,10))
plt.subplots_adjust(wspace=0.1, hspace=0.2)
axes = axes.flatten()
for image, title, axis in zip(images, titles, axes):
fig = sns.heatmap(np.flipud(image), cmap="Blues", ax=axis, xticklabels=range(1, 11),
yticklabels=list(reversed(range(12, 22))))
fig.set_ylabel('玩家点数总和', fontsize=15)
fig.set_xlabel('庄家明牌点数', fontsize=15)
fig.set_title(title, fontsize=15)
plt.show()
exit(0)
#print(V)
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(121)
sns.heatmap(np.flipud(V[12:,1:,0]), cmap="Wistia", xticklabels=range(1, 11),
yticklabels=list(reversed(range(12, 22))))
ax.set_ylabel('玩家的手牌', fontsize=16)
ax.set_xlabel('庄家的明牌', fontsize=16)
ax.set_title('无可用的 A')
ax = fig.add_subplot(122)
sns.heatmap(np.flipud(V[12:,1:,1]), cmap="Wistia", xticklabels=range(1, 11),
yticklabels=list(reversed(range(12, 22))))
ax.set_ylabel('玩家的手牌', fontsize=16)
ax.set_xlabel('庄家的明牌', fontsize=16)
ax.set_title('有可用的 A')
plt.show()
fig = plt.figure()
ax = Axes3D(fig)
ax.set_zticks(np.arange(-1,1,1))
X, Y = np.meshgrid(range(0,10), range(0,10))
ax.plot_surface(X, Y, V[12:,1:,0])
plt.show()

Просмотреть файл

@ -0,0 +1,182 @@
import numpy as np
import gym
import Algorithm.Base_MC_Policy_Iteration as base
import matplotlib.pyplot as plt
import seaborn as sns
import tqdm
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
HIT = 1
STICK = 0
def static_policy(state):
if state[0] >= 20:
return STICK
else:
return HIT
def run_V(env, gamma, episodes):
# 玩家手牌点数和1-21, 单元 0 空置
# 庄家明牌点数: 1-10, 单元 0 空置
# 有无可用的A 0-1 (无/有)
Value = np.zeros((22, 11, 2)) # G 的总和
Count = np.zeros((22, 11, 2)) # G 的数量
for episode in tqdm.trange(episodes):
s = env.reset()
Episode = [] # 一幕内的(状态,动作,奖励)序列
done = False
while (done is False): # 幕内循环
int_s = (s[0], s[1], int(s[2]))
action = static_policy(s)
next_s, reward, done, _ = env.step(action)
Episode.append((int_s, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, r = Episode[t]
G = gamma * G + r
Value[s] += G # 值累加
Count[s] += 1 # 数量加 1
print(Count)
print(np.min(Count[12:,1:,]))
Count[Count==0] = 1 # 把分母为0的填成1主要是针对终止状态Count为0
V = Value / Count # 求均值
return V
def greedy_policy(state, Q):
max_A = np.max(Q[state]) # 从几个动作中取最大值
argmax_A = np.where(Q[state] == max_A)[0]
action = np.random.choice(argmax_A)
return action
def run_Q_ES(env, gamma, episodes, policy):
# 玩家手牌点数和1-21, 单元 0 空置
# 庄家明牌点数: 1-10, 单元 0 空置
# 有无可用的A 0-1 (无/有)
# 动作0-停牌1-要牌
Value = np.zeros((22, 11, 2, 2)) # G 的总和
Count = np.ones((22, 11, 2, 2)) # G 的数量
for episode in tqdm.trange(episodes):
s = env.reset()
Episode = [] # 一幕内的(状态,动作,奖励)序列
# ES - initial random choice
int_s = (s[0], s[1], int(s[2]))
action = np.random.choice(2) # 这里必须随机取不能使用policy才符合探索的含义
next_s, reward, done, _ = env.step(action)
Episode.append((int_s, action, reward))
s = next_s # 迭代
while (done is False): # 幕内循环
int_s = (s[0], s[1], int(s[2])) # 把返回的 False/True 变成 0/1
action = np.random.choice(2, p=policy[int_s])
next_s, reward, done, _ = env.step(action)
Episode.append((int_s, action, reward))
s = next_s # 迭代
num_step = len(Episode)
G = 0
# 从后向前遍历计算 G 值
for t in range(num_step-1, -1, -1):
s, a, r = Episode[t]
G = gamma * G + r
Value[s][a] += G # 值累加
Count[s][a] += 1 # 数量加 1
# greedy policy
qs = Value[s] / Count[s]
policy[s] = 0
argmax = np.argmax(qs)
policy[s][argmax] = 1
return Value/Count
def create_blackjack_policy(args):
assert(args[0]+args[1]==1)
policy= np.zeros((22,11,2,2))
policy[:,:,:,0] = args[0]
policy[:,:,:,1] = args[1]
return policy
if __name__=="__main__":
policy = create_blackjack_policy((0.5,0.5))
gamma = 1
episodes = 500000
env = gym.make('Blackjack-v1', sab=True)
Q = run_Q_ES(env, gamma, episodes, policy)
env.close()
state_value_no_usable_ace = np.max(Q[12:, 1:, 0, :], axis=-1)
state_value_usable_ace = np.max(Q[12:, 1:, 1, :], axis=-1)
# get the optimal policy
action_no_usable_ace = np.argmax(Q[12:, 1:, 0, :], axis=-1)
action_usable_ace = np.argmax(Q[12:, 1:, 1, :], axis=-1)
print(Q[18,:,1,:])
images = [action_usable_ace,
state_value_usable_ace,
action_no_usable_ace,
state_value_no_usable_ace]
titles = ['最优策略(有可用的 A',
'最优状态价值(有可用的 A',
'最优策略(没有可用的 A',
'最优状态价值(没有可用的 A']
_, axes = plt.subplots(2, 2, figsize=(15,10))
plt.subplots_adjust(wspace=0.1, hspace=0.2)
axes = axes.flatten()
for image, title, axis in zip(images, titles, axes):
fig = sns.heatmap(np.flipud(image), cmap="Blues", ax=axis, xticklabels=range(1, 11),
yticklabels=list(reversed(range(12, 22))))
fig.set_ylabel('玩家点数总和', fontsize=15)
fig.set_xlabel('庄家明牌点数', fontsize=15)
fig.set_title(title, fontsize=15)
plt.show()
exit(0)
#print(V)
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(121)
sns.heatmap(np.flipud(V[12:,1:,0]), cmap="Wistia", xticklabels=range(1, 11),
yticklabels=list(reversed(range(12, 22))))
ax.set_ylabel('玩家的手牌', fontsize=16)
ax.set_xlabel('庄家的明牌', fontsize=16)
ax.set_title('无可用的 A')
ax = fig.add_subplot(122)
sns.heatmap(np.flipud(V[12:,1:,1]), cmap="Wistia", xticklabels=range(1, 11),
yticklabels=list(reversed(range(12, 22))))
ax.set_ylabel('玩家的手牌', fontsize=16)
ax.set_xlabel('庄家的明牌', fontsize=16)
ax.set_title('有可用的 A')
plt.show()
fig = plt.figure()
ax = Axes3D(fig)
ax.set_zticks(np.arange(-1,1,1))
X, Y = np.meshgrid(range(0,10), range(0,10))
ax.plot_surface(X, Y, V[12:,1:,0])
plt.show()

Просмотреть файл

@ -0,0 +1,93 @@
from enum import Enum
import numpy as np
class SeperatorLines(Enum):
empty = 0 # 打印空行
short = 1 # 打印10个'-'
middle = 2 # 打印30个'-'
long = 3 # 打印40个'='
def print_seperator_line(style: SeperatorLines, info=None):
if style == SeperatorLines.empty:
print("")
elif style == SeperatorLines.short:
if (info is None):
print("-"*10)
else:
print("----- " + info + " -----")
elif style == SeperatorLines.middle:
if (info is None):
print("-"*30)
else:
print("-"*15 + info + "-"*15)
elif style == SeperatorLines.long:
if (info is None):
print("="*40)
else:
print("="*20 + info + "="*20)
def print_V(dataModel, V):
vv = np.around(V,2)
print("状态价值函数计算结果(数组) :", vv)
for s in dataModel.S:
print(str.format("{0}:\t{1}", s.name, vv[s.value]))
def RMSE(x, y):
err = np.sqrt(np.sum(np.square(x - y))/y.shape[0])
return err
def test_policy(env, policy, episodes=100):
R = 0
for i in range(episodes):
s = env.reset()
done = False
while (done is False): # 幕内循环
action = np.argmax(policy[s])
# action = np.random.choice(env.action_space.n, p=policy[s])
next_s, reward, done, info = env.step(action)
R += reward
if done == True and reward == 0:
print(s, action, next_s)
s = next_s
return R
def create_policy(nS, nA, args):
assert(nA == len(args))
#shape = nS + (nA,)
shape = (nS, nA)
policy = np.zeros(shape)
sum = 0
for i in range(nA):
sum += args[i]
policy[:][i] = args[i]
assert(sum == 1)
return policy
# 从Q函数表格中抽取策略
def extract_policy_from_Q(Q, end_states):
policy = np.zeros_like(Q)
for s in range(Q.shape[0]):
if s not in end_states:
max_v = np.max(Q[s])
for a in range(Q[s].shape[0]):
if Q[s,a] == max_v:
policy[s, a] = 1
return policy
def calculat_V_from_Q(Q, policy):
nS = Q.shape[0]
V = np.zeros(Q.shape[0])
for s in range(nS):
V[s] = np.dot(policy[s], Q[s])
return V
if __name__=="__main__":
policy = create_policy(3, 2, (0.2,0.8))
print(policy)
print(policy[0,1])

Просмотреть файл

@ -0,0 +1,74 @@
import numpy as np
LEFT, DOWN, RIGHT, UP = 0, 1, 2, 3
LEFT_ARROW = 0x25c4
UP_ARROW = 0x25b2
RIGHT_ARROW = 0x25ba
DOWN_ARROW = 0x25bc
EMPTY_SPACE = 0x0020
CENTER_CROSS = 0x253c
SEP_LINE = 0x2500
class GridCell(object):
def __init__(self, q):
self.space = np.zeros((3,5), dtype=int)
self.space.fill(EMPTY_SPACE) # 空格
self.space[1,2] = CENTER_CROSS # 中心
self.q = np.round(q, 4)
# policy 是一个1x4的数组或矩阵,如[0.1,0.1,0.4,0.4]
# 这种情况要在向上和向右的地方画两个箭头
# 01234
# 0 ^
# 1 o >
# 2
if np.sum(q) != 0:
best_actions = np.argwhere(self.q == np.max(self.q))
pos = best_actions.flatten().tolist()
for action in pos:
if action == LEFT: # left
self.space[1,0] = LEFT_ARROW
self.space[1,1] = SEP_LINE
elif action == UP: # up
self.space[0,2] = UP_ARROW
elif action == RIGHT: # right
self.space[1,3] = SEP_LINE
self.space[1,4] = RIGHT_ARROW
elif action == DOWN: # down
self.space[2,2] = DOWN_ARROW
class Grid(object):
def __init__(self, Q, shape):
self.array = np.zeros((shape[0]*3, shape[1]*5), dtype=int)
for i in range(len(Q)):
row = (int)(i / shape[1])
col = (int)(i % shape[1])
q = Q[i]
cell = GridCell(q)
self.array[row*3:row*3+3, col*5:col*5+5] = cell.space
def draw(Q, shape):
grid = Grid(Q, shape)
for j, rows in enumerate(grid.array):
if (j % 3 == 0):
print("+-----"*shape[1], end="")
print("+")
print("|", end="")
for i,col in enumerate(rows):
print(chr(col), end="")
if ((i+1) % 5 == 0):
print("|", end="")
print()
print("+-----"*shape[1])
if __name__=="__main__":
Q = np.array([
[0.0155, 0.0164, 0.012, 0.015],
[0.0, 0.0, 0.00, 0.00],
[5.07, 3.06, 7.86 , 2.07],
[8.73, 8.73, 8.73 , 8.73]
])
draw(Q, (2,2))

Просмотреть файл

@ -0,0 +1,137 @@
import numpy as np
import matplotlib.pyplot as plt
from gym import spaces
LEFT, DOWN, RIGHT, UP = 0, 1, 2, 3
class GridWorld(object):
# 生成环境
def __init__(self,
GridWidth, GridHeight, StartStates, EndStates,
Actions, Policy, Transition,
StepReward, SpecialReward,
SpecialMove, Blocks):
self.Width = GridWidth
self.Height = GridHeight
self.Actions = Actions
self.nS = GridHeight * GridWidth
self.nA = len(Actions)
self.action_space = spaces.Discrete(self.nA)
self.observation_space = spaces.Discrete(self.nS)
self.SpecialReward = SpecialReward
self.StartStates = StartStates
self.EndStates = EndStates
self.SpecialMove = SpecialMove
self.Blocks = Blocks
self.Policy = self.__init_policy(Policy)
self.P = self.__init_states(Transition, StepReward)
# 把统一的policy设置复制到每个状态上
def __init_policy(self, Policy):
PI = {}
for s in range(self.nS):
PI[s] = Policy
return PI
# 用于生成状态->动作->转移->奖励字典
def __init_states(self, Transition, StepReward):
P_S_R_D = {}
s_id = 0
self.Pos2Sid = {}
self.Sid2Pos = {}
for y in range(self.Height):
for x in range(self.Width):
self.Pos2Sid[x,y] = s_id
self.Sid2Pos[s_id] = [x,y]
s_id += 1
for s, (x,y) in self.Sid2Pos.items():
P_S_R_D[s] = {}
if (s in self.EndStates):
continue
for action in self.Actions:
list_probs = []
for dir, prob in enumerate(Transition):
if (prob == 0.0):
continue
s_next = self.__get_next_state(
s, x, y, action + dir - 1) # 处理每一个转移概率方向逆时针减1
reward = StepReward # 通用奖励定义 (-1)
if (s, s_next) in self.SpecialReward: # 如果有特殊奖励定义
reward = self.SpecialReward[(s, s_next)]
list_probs.append((prob, s_next, reward, self.is_end(s_next)))
P_S_R_D[s][action] = list_probs
return P_S_R_D
# 用于计算移动后的下一个状态
# 左上角为 [0,0], 横向为 x, 纵向为 y
def __get_next_state(self, s, x, y, action):
action = action % 4 # 避免负数
if (s,action) in self.SpecialMove:
return self.SpecialMove[(s,action)]
if (action == UP): # 向上转移
if (y != 0): # 不在上方边界处,否则停在原地不动
s = s - self.Width
elif (action == DOWN): # 向下转移
if (y != self.Height-1):# 不在下方边界处,否则停在原地不动
s = s + self.Width
elif (action == LEFT): # 向左转移
if (x != 0): # 不在左侧边界处,否则停在原地不动
s = s - 1
elif (action == RIGHT): # 向右转移
if (x != self.Width-1): # 不在右侧边界处,否则停在原地不动
s = s + 1
return s
def is_end(self, s):
return (s in self.EndStates)
def reset(self):
self.curr_state = np.random.choice(self.StartStates)
return self.curr_state
def close(self):
pass
def step(self, a):
transitions = self.P[self.curr_state][a]
num = len(transitions)
if (num == 1):
self.curr_state = transitions[0][1]
return transitions[0][1], transitions[0][2], transitions[0][3], transitions[0][0]
else:
self.p = [transitions[i][0] for i in range(num)]
item = np.random.choice(self.nA, p=self.p)
self.curr_state = item[0][1]
return item[1], item[2], item[3], item[0]
action_names = ['LEFT', 'UP', 'RIGHT', 'DOWN']
def print_P(P):
print("状态->动作->转移->奖励 字典:")
for s,v in P.items():
print("state =",s)
for action,v2 in v.items():
print(str.format("\taction = {0} : {1}", action_names[action], v2))
# left, up, right, down
chars = [0x2190, 0x2191, 0x2192, 0x2193]
# 需要处理多个值相等的情况
def print_policy(policy, shape):
best_actions = np.argmax(policy, axis=1)
for i, action in enumerate(best_actions):
print(chr(chars[action]), end="")
print(" ", end="")
if ((i+1) % shape[0] == 0):
print("")
# 绘图
def draw_table(V, shape):
tab = plt.table(cellText=V, loc='center', rowHeights=[0.1]*5)
tab.scale(1,1)
plt.axis('off')
plt.show()

Двоичный файл не отображается.

Просмотреть файл

@ -14,9 +14,9 @@
- 为什么叫单臂/多臂强盗?
- 因为这机器回报率低,但是容易上瘾,让玩家输个精光,所以叫一条胳膊的强盗。很多现代的老虎机仍然保持一个传统的操作摇臂,可与按钮并行使用。
- 因为这机器回报率低,但是容易上瘾,让玩家输个精光,所以叫一条胳膊的强盗。很多现代的老虎机仍然保持一个传统的机械操作摇臂,可与电子按钮并行使用。
- 有的赌徒有策略的操作多台老虎机,以实现最高的可能获利这样的赌徒就称作“多臂强盗”K-Armed Bandits
- 有的赌徒有策略地操作多台老虎机,以实现最高的获利这样的赌徒就称作“多臂强盗”K-Armed Bandits
- 某些翻译软件把 K-Armed Bandits 翻译成“武装匪徒”,因为 Armed 有“武装”的意思,读者自行理解吧。
@ -26,21 +26,21 @@
1. Choose Slot Machines with High Price per Spin
选择每次摇臂操作赌注高的老虎机。
选择每次赌注高的老虎机。
虽然这一提示似乎违反直觉,但下注的赌注越高的老虎机回报的频率越高,金额越大。这就是老虎机的设计方式,吐币的百分比与赌注大小成正比。此外,最大赌注可以帮助您解锁特殊功能,如免费旋转或头奖
虽然这一提示似乎违反直觉,但下注的赌注越高的老虎机回报的频率越高,金额越大。这就是老虎机的设计方式,吐币的百分比与赌注大小成正比。此外,最大赌注可以帮助您解锁特殊功能,如免费增加一次摇臂机会或中头奖机会
2. Choose the Least Complicated Slot Machines
选择最简单的老虎机。
在赌场使用这条普遍的经验法则:如果游戏复杂,它的回报就会少。在复杂的老虎机游戏中,赌场试图通过特殊的乘数、奖金或额外的生命来吸引你的注意力。虽然您可能喜欢这些特殊功能,但它们只会使游戏复杂化,并降低获胜的几率。
在赌场使用这条普遍的经验法则:如果游戏复杂,它的回报就会少。在复杂的老虎机游戏中,赌场试图通过特殊的乘数、奖金或额外的摇臂次数来吸引你的注意力。虽然你可能喜欢这些特殊功能,但它们只会使游戏复杂化,并降低获胜的几率。
3. Test the Machines Before you Play
在玩之前测试机器。
在玩大赌注之前测试机器。
大多数赌场在其许多老虎机上提供免费试用,这意味着你可以先测试五到十圈,让你更容易制定出一个行动计划,以便在使用真金白银的时候采取行动
大多数赌场在其许多老虎机上提供免费试用,这意味着你可以先测试五到十圈,让你更容易制定出一个行动计划,以便在使用真金白银的时候采取正确的策略
4. Use Cash Instead of Card
@ -67,8 +67,8 @@
因为单臂的赌博机应用场景比较简单,所以在算法研究领域都是多臂赌博机平台,有以下几类:
- 随机多臂赌博机
- 贝叶斯多臂赌博机
- 固定收益的多臂赌博机
- 随机分布收益的多臂赌博机
- 上下文赌博机
- 其它
@ -78,7 +78,7 @@
#### 固定收益的多臂赌博机
如图 2.1.1 所示,假设是一个三臂赌博机,玩家投一次币(假设是一元钱),选择一个序号,拉动一次拉杆臂(或者按一次按钮)。
如图 2.1.1 所示,假设是一个三臂赌博机,玩家投一次币(假设是一元钱),选择一个序号,拉动一次臂(或者按一次按钮)。
<center>
@ -87,14 +87,14 @@
图 2.1.1 固定收益的多臂赌博机
</center>
假设玩家选择 1 号,则该臂机会以 $p_1=0.3$ 的概率回吐一元钱做为奖励,即 10 次中可能有 3 次回吐一元钱,其它 7 次没有奖励。另外两臂的回吐概率是 $p_2=0.5,p_3=0.8$,但是玩家在开始时并不知道这个信息,需要摸索。
假设玩家选择 1 号,则该臂机会以 $p_1=0.3$ 的概率回吐一元钱做为奖励,即 10 次中可能有 3 次回吐一元钱,其它 7 次没有奖励。另外两臂的回吐概率是 $p_2=0.5,p_3=0.8$,但是玩家在开始时**并不知道**这个信息,需要摸索。
这个应用场景相对简单,可以总结为:
- 收益固定为离散值 0/1即伯努利分布。
- 收益概率可以灵活设置。
- 回吐概率可以灵活设置。
- 算法简单,运行速度快,见效快。
- 初者容易理解。
- 初者容易理解。
- 有预设的上限可以参考:比如以图 2.1.1 的设置为例,玩 1000 次的最大收益应该是 800 元。
- 衡量指标有限,不能体现复杂算法的优势。
- 能体现客观世界(赌场)的真实情况,但是对算法研究没有帮助,所以不能应用于其它场景。
@ -112,7 +112,7 @@
图 2.1.2 收益按概率分布的多臂赌博机
</center>
依然假设只有三臂,每个拉杆臂的收益都是一个正态分布,方差为 1均值有可能是 0,-1,+1。
依然假设只有三臂,每个拉杆臂的收益都是一个正态分布,方差为 1均值有可能分别是 0,-1,+1。
举例来说,假设用户选择 1 号臂,进行 200 轮游戏的话,会看到如图 2.1.3 的奖励分布。
@ -129,8 +129,8 @@
- 奖励(即收益)不是正整数,而是浮点数,有正有负;
- 奖励在 0 附近的次数最多,占据了 70% 左右;
- 最大值,有机会得到 +2 的奖励或更大,次数极少;
- 最小值,有机会得到 -2 的奖励或更小,次数极少。
- 最大值,有机会得到 +2 的奖励或更大,次数极少;
- 最小值,有机会得到 -2 的奖励或更小,次数极少。
2 号和 3 号臂就是 1 号臂的分布情况分别向左或右移动一个单位的正态分布。

Просмотреть файл

@ -1,6 +1,6 @@
## 2.2 多臂赌博机建模
### 2.2.1
### 2.2.1 收益随机分布的多臂赌博机
下面再描述一下我们将要研究的多臂赌博机平台的细节。
@ -11,16 +11,16 @@
<img src='./img/One-Arm.png'/>
图 2.2.1 多臂赌博机的一个臂的奖励情况
(左图:奖励及其分布次数;右图:用小提琴图表示的)
(左图:奖励及其分布次数;右图:用小提琴图表示的收益分布
</center>
- 每个臂的奖励都是一个正态分布,均值为 0方差为 1样本数量为 2000。
- 最小值和最大值都是采样的结果,基本是在 $\pm3$ 左右。
- 最小值和最大值都是随机采样的结果,基本是在 $\pm3$ 左右。
- 均值在 0 左右,偏差很小。
- 小提琴蓝色区域的宽度表示样本的密度。
- 上下四分位和置信区间在这个问题中暂时用不上。
生成原始数据的代码如下:
生成 10 臂赌博机的原始数据的代码如下:
【代码位置】bandit_22_ArmBandits.py
@ -29,12 +29,12 @@
num_arm = 10
num_data = 2000
np.random.seed(5)
k_reward_dist = np.random.randn(num_data, num_arm) # num_arm 放在后面是为了可以做加法
k_reward_dist = np.random.randn(num_data, num_arm) # num_arm 放在后面是为了可以做矩阵加法
print("原始均值=", np.round(np.mean(k_reward_dist, axis=0),3))
draw_one_arm(k_reward_dist[:,0])
```
得到原始数据的均值为:
得到 10 臂的原始数据的均值为:
```
原始均值= [ 0.018 0.019 0.004 -0.02 -0.005 0.004 -0.011 -0.005 0.027 -0.015]
@ -108,10 +108,10 @@ $$
<img src='./img/K-arm-bandits.png'/>
图 2.2.3 多臂赌博机奖励分布
上图10臂赌博机期望数据分布下图10臂赌博机按均值排序分布
上图10 臂赌博机期望数据分布下图10 臂赌博机按均值排序)
</center>
图 2.2.3 的上图,是按期望均值的原始顺序排列的 10 臂赌博机的奖励分布情况,下图是排序后的分布情况,比如:上图中的 1 号臂,在下图中排在 5 号位置这样一来10 个臂的回报情况是从左到右依次变好。当然,我们在做算法的时候,要假装不知道 10 号臂是最佳选择,而是盲猜
图 2.2.3 的上图,是按期望均值的原始顺序排列的 10 臂赌博机的奖励分布情况,下图是排序后的分布情况,比如:上图中的 1 号臂,在下图中排在 5 号位置这样一来10 个臂的回报情况是从左到右依次变好。当然,我们在做算法的时候,要假装不知道 10 号臂是最佳选择,而是需要探索
### 2.2.2 如何定义最好的动作
@ -137,7 +137,7 @@ q(a_3)&=\frac{1.3+1.0}{2}=1.15
\tag{2.2.1}
$$
2.2.1)告诉我们,在只玩 10 轮的经验下,$q(a_1)$ 的值最大,$a_1$ 就是最好的动作。
2.2.1)告诉我们,在只玩 10 轮的经验下,$q(a_3)$ 的值最大,$a_3$ 就是最好的动作。
理论上,可以定义某个动作的价值为 $q_*$,则有:
@ -158,7 +158,7 @@ $n$ 表示动作被选择的次数,当 $n$ 足够大时,可以保证每个
在式2.2.1)中,当前状态对于 $a_3$ 来说是 $n=2$。假设 $t=11$ 时选择了 $a_3$,则 $n=3$$q(a_1),q(a_2)$ 的值不会变化,$q(a_3)$ 的值会被重新计算,由$q_2(a_3)$ 变成 $q_3(a_3)$。介绍一个小技巧,计算 $q_3(a_3)$ 时可以这样做:
$$
q_{3}(a_3)= \frac{1}{3}[2q_{2}(a_3) + R_{3}]=\frac{1}{3}[3q_{2}(a_3)-q_{2}+R_{3}]=q_{2}(a_3)+\frac{1}{3}[R_{3}-q_{2}(a_3)]
q_{3}(a_3)= \frac{1}{3}[2q_{2}(a_3) + R_{3}]=\frac{1}{3}[3q_{2}(a_3)-q_{2}(a_3)+R_{3}]=q_{2}(a_3)+\frac{1}{3}[R_{3}-q_{2}(a_3)]
\tag{2.2.4}
$$
@ -182,6 +182,6 @@ Q_{n}&=Q_{n-1}+\alpha [R_n - Q_{n-1}]
\tag{2.2.6}
$$
步长参数 $\alpha \in (0,1]$ 是一个常数而不是像式2.2.5)的 $\frac{1}{n}$ 那样随着采样次数的增加而降低。在解决一个非平稳问题的时候,即无法用 $\frac{1}{n}$ 来做平均值计算的时候可以使用式2.2.6)。
步长参数 $\alpha \in (0,1]$ 是一个常数而不是像式2.2.5)的 $\frac{1}{n}$ 那样随着采样次数 $n$ 的增加而降低。在解决一个非平稳问题的时候,即无法用 $\frac{1}{n}$ 来做平均值计算的时候可以使用式2.2.6)。
在 2.4 节贪心法中会使用式2.2.5);在 2.5 节梯度上升法中会使用式2.2.6)。

Некоторые файлы не были показаны из-за слишком большого количества измененных файлов Показать больше