蒙特卡洛方法（Monte-Carlo methods）也被称为统计模拟方法，是一种基于概率统计的数值计算方法。运用蒙特卡洛方法时，我们通常使用重复随机抽样，然后运用概率统计方法来从抽样结果中归纳出我们想求的目标的数值估计。
一个简单的例子是用蒙特卡洛方法来计算圆的面积。例如，在下图所示的正方形内部随机产生若干个点，细数落在圆中点的个数，圆的面积与正方形面积之比就等于圆中点的个数与正方形中点的个数之比。如果我们随机产生的点的个数越多，计算得到圆的面积就越接近于真实的圆的面积。
强化学习算法的精髓之一是解决无模型的马尔可夫决策问题。无模型的强化学习算法主要包括蒙特卡洛方法和时序差分方法。本文的主题是蒙特卡洛方法，那么如何用蒙特卡洛方法来估计一个策略在一个马尔可夫决策过程中的状态价值函数呢？
回到值函数最原始的定义公式
V
(
s
)
=
E
[
G
t
∣
S
t
=
s
]
V(s)=\mathbb{E}\left[G_{t} \mid S_{t}=s\right]
V(s)=E[Gt∣St=s]
状态值函数的计算实际上是计算返回值的期望。动态规划的方法是利用模型计算该期望。在没有模型时，我们可以采用蒙特卡洛的方法计算该期望，即利用随机样本估计期望。
一个很直观的想法就是用策略在 MDP 上采样很多条序列，计算从这个状态出发的回报，再求其期望就可以了，公式如下：
V
π
(
s
)
=
E
π
[
G
t
∣
S
t
=
s
]
≈
1
N
∑
i
=
1
N
G
t
(
i
)
V^{\pi}(s)=\mathbb{E}_{\pi}\left[G_{t} \mid S_{t}=s\right] \approx \frac{1}{N} \sum_{i=1}^{N} G_{t}^{(i)}
Vπ(s)=Eπ[Gt∣St=s]≈N1i=1∑NGt(i)
在计算值函数时，蒙特卡洛方法是利用经验平均代替随机变量的期望。下面分别阐述经验和平均的意思。
“经验”：当要评估智能体的当前策略时，我们可以利用策略产生很多次试验，每次试验都是从任意的初始状态开始直到终止，比如一次试验（episode）为
S
1
,
A
1
,
R
2
,
⋯
,
S
T
S_1,A_1,R_2,\cdots,S_T
S1,A1,R2,⋯,ST，计算一次试验中状态
s
s
s处的折扣回报返回值为
G
(
s
)
=
R
t
+
1
+
γ
R
t
+
2
+
⋯
+
γ
T
−
1
R
T
G(s)=R_{t+1}+\gamma R_{t+2}+\cdots+\gamma^{T-1}R_T
G(s)=Rt+1+γRt+2+⋯+γT−1RT，那么"经验"就是指利用该策略做很多次试验，产生很多幕数据，见下图。
“平均”：就是求平均值。不过，利用蒙特卡洛方法求状态处的值函数时，又可以分为第一次访问蒙特卡罗方法和每次访问蒙特卡罗方法。第一次访问蒙特卡洛方法是指在计算状态处的值函数时，只利用每次试验中第一次访问到状态 s s s时的返回值。如上图中第一次试验所示，计算状态 s s s处的均值时只利用 G 11 G_{11} G11，因此第一次访问蒙特卡洛方法的计算公式为
V ( s ) = G 11 ( s ) + G 21 ( s ) + ⋯ N ( s ) V(s)=\dfrac{G_{11}(s)+G_{21}(s)+\cdots}{N(s)} V(s)=N(s)G11(s)+G21(s)+⋯
每次访问蒙特卡洛方法是指在计算状态 s s s处的值函数时，利用所有访问到状态 s s s时的回报返回值，即
V ( s ) = G 11 ( s ) + G 12 ( s ) + ⋯ + G 21 ( s ) + ⋯ N ( s ) V(s)=\dfrac{G_{11}(s)+G_{12}(s)+\cdots+G_{21}(s)+\cdots}{N(s)} V(s)=N(s)G11(s)+G12(s)+⋯+G21(s)+⋯
同时，根据大数定律：当 N ( s ) → ∞ N(s)\rightarrow \infty N(s)→∞，有 V ( s ) → V π ( s ) V(s)\rightarrow V^{\pi}(s) V(s)→Vπ(s)。
在一条序列中，可能没有出现过这个状态，可能只出现过一次这个状态，也可能出现过很多次这个状态。我们介绍的蒙特卡洛价值估计方法会在该状态每一次出现时计算它的回报。还有一种选择是一条序列只计算一次回报，也就是这条序列第一次出现该状态时计算后面的累积奖励，而后面再次出现该状态时，该状态就被忽略了。假设我们现在用策略从状态开始采样序列，据此来计算状态价值。我们为每一个状态维护一个计数器和总回报，计算状态价值的具体过程如下所示。
由于智能体与环境交互的模型是未知的，蒙特卡罗方法是利用经验平均来估计值函数，而能否得到正确的值函数，则取决于经验。因此，如何获得充足的经验是无模型强化学习的核心所在。
在动态规划方法中，为了保证值函数的收敛性，算法会逐个扫描状态空间中的状态。而无模型的方法充分评估策略值函数的前提，是每个状态都能被访问到。因此，在蒙特卡洛方法中必须采用一定的方法保证每个状态都能被访问到，方法之一是探索性初始化。探索性初始化是指每个状态都有一定的几率作为初始状态。在学习基于探索性初始化的蒙特卡罗方法前，我们还需要先了解策略改善方法，以及便于进行迭代计算的平均方法。
蒙特卡洛策略改善方法：蒙特卡洛方法利用经验平均估计策略值函数。估计出值函数后，对于每个状态，它通过最大化动作值函数来进行策略的改善。即
π
(
s
)
=
arg
max
a
q
(
s
,
a
)
\pi(s)=\arg\max_a q(s,a)
π(s)=argamaxq(s,a)
增量更新方法：计算回报的期望时，除了可以把所有的回报加起来除以次数，还有一种增量更新的方法。对于每个状态和对应回报，进行如下计算：
探索性初始化在迭代每一幕时，初始状态是随机分配的，这样可以保证迭代过程中每个状态行为对都能被选中。它蕴含着一个假设：假设所有的动作都被无限频繁选中。对于这个假设，有时很难成立，或无法完全保证。所以，需要设计一个温和的探索策略，即对所有的状态 s s s和 a a a都满足： π ( a ∣ s ) > 0 \pi(a|s)>0 π(a∣s)>0。也就是说，温和的探索策略是指在任意状态下，采用动作集中每个动作的概率都大于零。典型的温和策略是 ε \varepsilon ε-soft策略：
π ( a ∣ s ) ← { 1 − ε + ε ∣ A ( s ) ∣ if a = arg max a Q ( s , a ) ε ∣ A ( s ) ∣ if a ≠ arg max a Q ( s , a ) \pi(a \mid s) \leftarrow\left\{\begin{array}{ll} 1-\varepsilon+\frac{\varepsilon}{|A(s)|} & \text { if } a=\arg \max _{a} Q(s, a) \\ \frac{\varepsilon}{|A(s)|} & \text { if } a \neq \arg \max _{a} Q(s, a) \end{array}\right. π(a∣s)←{1−ε+∣A(s)∣ε∣A(s)∣ε if a=argmaxaQ(s,a) if a=argmaxaQ(s,a)
根据探索策略（行动策略）和评估的策略是否为同一个策略，蒙特卡罗方法又分为on-policy和off-policy两种方法。若行动策略和评估及改善的策略是同一个策略，称为on-policy，可翻译为同策略。比如，要产生数据的策略和评估及要改善的策略都是 ε \varepsilon ε-soft策略。若行动策略和评估及改善的策略是不同的策略，称为off-policy，可翻译为异策略。异策略可以保证充分的探索性。例如用来评估和改善的策略 π \pi π是贪婪策略，用于产生数据的探索性策略 μ \mu μ为探索性策略，如 ε \varepsilon ε-soft策略。
用于异策略的目标策略 π \pi π和行动策略 μ \mu μ并非任意选择的，而是必须满足一定的条件。这个条件是覆盖性条件，即行动策略 μ \mu μ产生的行为覆盖或包含目标策略 π \pi π产生的行为。利用式子表示：满足 π ( a ∣ s ) > 0 \pi(a|s)>0 π(a∣s)>0的任何 ( s , a ) (s,a) (s,a)均满足 μ ( a ∣ s ) > 0 \mu(a|s)>0 μ(a∣s)>0。利用行为策略产生的数据评估目标策略需要利用重要性采样方法。
重要性采样来源于求期望，重要性采样的原理可以用图来描述。
E
[
f
]
=
∫
f
(
z
)
p
(
z
)
d
z
E[f]=\int f(z) p(z) dz
E[f]=∫f(z)p(z)dz
当随机变量
z
z
z的分布非常复杂时，无法利用解析的方法产生用于逼近期望的样本，这时，我们可以选用一个概率分布很简单，很容易产生样本的概率分布
q
(
z
)
q(z)
q(z)，比如正态分布。原来的期望可变为
E
[
f
]
=
∫
f
(
z
)
p
(
z
)
d
z
=
∫
f
(
z
)
p
(
z
)
q
(
z
)
q
(
z
)
d
z
=
1
N
∑
n
p
(
z
n
)
q
(
z
n
)
f
(
z
n
)
,
z
n
∼
q
(
z
)
\begin{aligned} E[f] &=\int f(z) p(z) d z \\ &=\int f(z) \frac{p(z)}{q(z)} q(z) d z \\ &=\frac{1}{N} \sum_{n} \frac{p\left(z^{n}\right)}{q\left(z^{n}\right)} f\left(z^{n}\right), z^{n} \sim q(z) \end{aligned}
E[f]=∫f(z)p(z)dz=∫f(z)q(z)p(z)q(z)dz=N1n∑q(zn)p(zn)f(zn),zn∼q(z)
定义重要性权重：
w
n
=
p
(
z
n
)
q
(
z
n
)
w^n=\dfrac{p(z^n)}{q(z^n)}
wn=q(zn)p(zn)，普通的重要性采样求积分如下式所示：
E
[
f
]
=
1
N
∑
n
ω
n
f
(
z
n
)
E[f]=\frac{1}{N} \sum_{n} \omega^{n} f\left(z^{n}\right)
E[f]=N1n∑ωnf(zn)
通过上面的式子，可以分析出基于重要性采样的积分估计为无偏估计，即估计的期望值等于真实的期望。但是，基于重要性采样的积分估计的方差无穷大。这是因为原来的被积函数乘了一个重要性权重，改变了被积函数的形状及分布。尽管被积函数的均值没有发生变化，但方差明显发生改变。
在重要性采样中，使用的采样概率分布与原概率分布越接近，方差越小。然而，被积函数的概率分布往往很难求得，因此没有与之相似的简单采样概率分布，如果使用分布差别很大的采样概率对原概率分布进行采样，方差会趋近于无穷大。一种减小重要性采样积分方差的方法是采用加权重要性采样：
E
[
f
]
≈
∑
n
=
1
N
ω
n
∑
m
=
1
N
ω
m
f
(
z
n
)
E[f] \approx \sum_{n=1}^{N} \dfrac{\omega^{n}}{\sum_{m=1}^{N} \omega^{m}} f\left(z^{n}\right)
E[f]≈n=1∑N∑m=1Nωmωnf(zn)
在异策略方法中，行动策略 μ \mu μ即用来产生样本的策略，所产生的轨迹概率分布相当于重要性采样中的 q ( z ) q(z) q(z)， 用来评估和改进的策略 π \pi π所对应的轨迹概率分布为 p ( z ) p(z) p(z)，因此利用行动策略 μ \mu μ所产生的累积函数返回值来评估策略 π \pi π时，需要在累积函数返回值前面乘以重要性权重。
在目标策略
π
\pi
π下，一次试验的概率为
P
r
(
A
t
,
S
t
+
1
,
⋯
,
S
T
)
=
∏
k
=
t
T
−
1
π
(
A
k
∣
S
k
)
p
(
S
k
+
1
∣
S
K
,
A
k
)
Pr(A_t,S_{t+1},\cdots,S_T)=\prod_{k=t}^{T-1}\pi(A_k|S_k)p(S_{k+1}|S_K,A_k)
Pr(At,St+1,⋯,ST)=k=t∏T−1π(Ak∣Sk)p(Sk+1∣SK,Ak)
在行动策略
μ
\mu
μ下，相应的试验的概率为
P
r
(
A
t
,
S
t
+
1
,
⋯
,
S
T
)
=
∏
k
=
t
T
−
1
μ
(
A
k
∣
S
k
)
p
(
S
k
+
1
∣
S
K
,
A
k
)
Pr(A_t,S_{t+1},\cdots,S_T)=\prod_{k=t}^{T-1}\mu(A_k|S_k)p(S_{k+1}|S_K,A_k)
Pr(At,St+1,⋯,ST)=k=t∏T−1μ(Ak∣Sk)p(Sk+1∣SK,Ak)
因此重要性权重为
ρ
t
T
=
∏
k
=
t
T
−
1
π
(
A
k
∣
S
k
)
p
(
S
k
+
1
∣
S
k
,
A
k
)
∏
k
=
1
T
−
1
μ
(
A
k
∣
S
k
)
p
(
S
k
+
1
∣
S
k
,
A
k
)
=
∏
k
=
t
T
−
1
π
(
A
k
∣
S
k
)
μ
(
A
k
∣
S
k
)
\rho_{t}^{T}=\frac{\prod_{k=t}^{T-1} \pi\left(A_{k} \mid S_{k}\right) p\left(S_{k+1} \mid S_{k}, A_{k}\right)}{\prod_{k=1}^{T-1} \mu\left(A_{k} \mid S_{k}\right) p\left(S_{k+1} \mid S_{k}, A_{k}\right)}=\prod_{k=t}^{T-1} \frac{\pi\left(A_{k} \mid S_{k}\right)}{\mu\left(A_{k} \mid S_{k}\right)}
ρtT=∏k=1T−1μ(Ak∣Sk)p(Sk+1∣Sk,Ak)∏k=tT−1π(Ak∣Sk)p(Sk+1∣Sk,Ak)=k=t∏T−1μ(Ak∣Sk)π(Ak∣Sk)
普通重要性采样的值函数估计式如下
V ( s ) = ∑ t ∈ T ( s ) ρ t T ( t ) G t ∣ T ( s ∣ V(s)=\dfrac{\sum_{t\in \mathcal{T}(s)}\rho_t^{T(t)}G_t}{|\mathcal{T}(s|} V(s)=∣T(s∣∑t∈T(s)ρtT(t)Gt
加权重要性采样值函数估计为
V ( s ) = ∑ t ∈ T ( s ) ρ t T ( t ) G t ∑ t ∈ T ( s ) ρ t T ( t ) V(s)=\dfrac{\sum_{t\in \mathcal{T}(s)}\rho_t^{T(t)}G_t}{\sum_{t\in \mathcal{T}(s)}\rho_t^{T(t)}} V(s)=∑t∈T(s)ρtT(t)∑t∈T(s)ρtT(t)Gt
根据上述分析，异策略每次访问蒙特卡罗算法的伪代码见如下算法伪代码。
Input：对于所有的
s
∈
S
s \in S
s∈S,
a
∈
A
(
s
)
a \in A(s)
a∈A(s)：
(
s
,
a
)
←
(s,a) ←
(s,a)←任意，
C
(
s
,
a
)
←
0
C(s,a) ← 0
C(s,a)←0,
π
(
s
)
←
π(s) ←
π(s)←相对于
Q
Q
Q的贪婪策略
Repeat：
Until
与基于动态规划的方法相比，基于蒙特卡罗的方法只是在值函数估计上有所不同，在整个框架上则是相同的，即估当前策略，再利用学到的值函数进行策略改善。
# 把输入的两个字符串通过“-”连接
def join(str1, str2):
return str1 + '-' + str2
def sample(MDP, Pi, timestep_max, number):
''' 采样函数,策略Pi,限制最长时间步timestep_max,总共采样序列数number '''
S, A, P, R, gamma = MDP
episodes = []
for _ in range(number):
episode = []
timestep = 0
s = S[np.random.randint(4)] # 随机选择一个除s5以外的状态s作为起点
# 当前状态为终止状态或者时间步太长时,一次采样结束
while s != "s5" and timestep <= timestep_max:
timestep += 1
rand, temp = np.random.rand(), 0
# 在状态s下根据策略选择动作
for a_opt in A:
temp += Pi.get(join(s, a_opt), 0)
if temp > rand:
a = a_opt
r = R.get(join(s, a), 0)
break
rand, temp = np.random.rand(), 0
# 根据状态转移概率得到下一个状态s_next
for s_opt in S:
temp += P.get(join(join(s, a), s_opt), 0)
if temp > rand:
s_next = s_opt
break
episode.append((s, a, r, s_next)) # 把（s,a,r,s_next）元组放入序列中
s = s_next # s_next变成当前状态,开始接下来的循环
episodes.append(episode)
return episodes
# 采样5次,每个序列最长不超过1000步
episodes = sample(MDP, Pi_1, 20, 5)
print('第一条序列\n', episodes[0])
print('第二条序列\n', episodes[1])
print('第五条序列\n', episodes[4])
# 对所有采样序列计算所有状态的价值
def MC(episodes, V, N, gamma):
for episode in episodes:
G = 0
for i in range(len(episode) - 1, -1, -1): #一个序列从后往前计算
(s, a, r, s_next) = episode[i]
G = r + gamma * G
N[s] = N[s] + 1
V[s] = V[s] + (G - V[s]) / N[s]
timestep_max = 20
# 采样1000次,可以自行修改
episodes = sample(MDP, Pi_1, timestep_max, 1000)
gamma = 0.5
V = {"s1": 0, "s2": 0, "s3": 0, "s4": 0, "s5": 0}
N = {"s1": 0, "s2": 0, "s3": 0, "s4": 0, "s5": 0}
MC(episodes, V, N, gamma)
print("使用蒙特卡洛方法计算MDP的状态价值为\n", V)
\quad
\quad
参考
持续更新中…