蒙特卡洛方法_蒙特卡罗求圆周率-CSDN博客

本文链接：https://blog.csdn.net/lesileqin/article/details/121534911

文章目录

蒙特卡洛方法（Monte-Carlo methods）
- 1）计算 $\pi$
- 2）计算状态价值函数
参考文献

蒙特卡洛方法（Monte-Carlo methods）

1）计算 $\pi$

蒙特卡洛方法也被称为统计模拟方法，是一种基于概率统计的数值计算方法。

下面以计算 $\pi$ 为例，介绍蒙特卡洛方法：
在这里插入图片描述

例如在上图的一个边长为2的正方形中又一个圆（暂且不管里面的点），我们可以很轻松的得到：

正方形的面积 $A_1 = 2^2 = 4$
圆的面积 $A_2 = \pi r^2 = \pi 1^2 = \pi$

现在有无数个点均匀分布在正方形内，如上图所示，那么可以轻松得到点落在圆内的概率是： $\frac{A_2}{A_1} = \frac{\pi}{4}$

在正方形内随机抽样n个点，所以点落在圆内的期望 $P_n=\frac{\pi n}{4}$

假设有m个点落在圆内，如果m非常大的话，就有： $\approx \frac{n \pi}{4}$ ，即 $\pi \approx \frac{4m}{n}$

于是就有了算法：

User specify a big n ; reset counter m = 0
For i = 1 to n:
	a) Randomly generate x belong [-1,1]
	b) Randomly generate y belong [-1,1]
	c) if x*x + y*y <= 1 , then m = m + 1
Return pi approx (4m) / n

使用python代码实现：

import numpy as np

n = 10000000
m = 0
for i in range(n):
    x = (1-(-1)) * np.random.random() + (-1)
    y = (1-(-1)) * np.random.random() + (-1)
    if x*x + y*y <= 1:
        m = m + 1
(4 * m) / n

-------------------------------------------------------
3.1424084

结果还不错，如果n无限大，就会越来越接近 $\pi$

2）计算状态价值函数

使用蒙特卡洛方法来估计一个策略在一个马尔科夫决策过程中的状态价值函数。

用策略在MDP上采样很多条序列，然后计算从这个状态出发的回报的再求平均值：
在这里插入图片描述
步骤：

使用策略 $\pi$ 采样：
对每一条序列中的每一时间步 t 的状态 s 进行以下操作：
- 更新状态 s 的计数器： $N (s) = N (s) + 1$
- 更新状态 s 的总回报： $M(s) = M(s) + G_t$
为每一个状态的价值被估计为回报的平均值： $V (s) = M (s) / N (s)$

还有一种增量更新的方法：对于每个状态 s 和对应回报 G ，进行：

$N (s) = N (s) + 1$
$\frac{1}{N(s)}(G-V(s))$

一个小例子：
在这里插入图片描述

使用python代码实现：

S = ["S1","S2","S3","S4","S5"] # 状态集合
A = ["保持S1","前往S1","前往S2","前往S3","前往S4","前往S5","概率前往"] # 动作集合
# 状态转移函数
P = {
    "S1-保持S1-S1":1.0, "S1-前往S2-S2":1.0,
    "S2-前往S1-S1":1.0, "S2-前往S3-S3":1.0,
    "S3-前往S4-S4":1.0, "S3-前往S5-S5":1.0,
    "S4-前往S5-S5":1.0, "S4-概率前往-S2":0.2,
    "S4-概率前往-S3":0.4, "S4-概率前往-S4":0.4,
}
# 奖励函数
R = {
    "S1-保持S1":-1, "S1-前往S2":0,
    "S2-前往S1":-1, "S2-前往S3":-2,
    "S3-前往S4":-2, "S3-前往S5":0,
    "S4-前往S5":10, "S4-概率前往":1,
}
gamma = 0.5 # 折扣因子
MDP = (S, A, P, R, gamma)

# 策略1，随机策略
Pi_1 = {
    "S1-保持S1":0.5, "S1-前往S2":0.5,
    "S2-前往S1":0.5, "S2-前往S3":0.5,
    "S3-前往S4":0.5, "S3-前往S5":0.5,
    "S4-前往S5":0.5, "S4-概率前往":0.5,
}
# 策略2
Pi_2 = {
    "S1-保持S1":0.6, "S1-前往S2":0.4,
    "S2-前往S1":0.3, "S2-前往S3":0.7,
    "S3-前往S4":0.5, "S3-前往S5":0.5,
    "S4-前往S5":0.1, "S4-概率前往":0.9,
}
# 把输入的两个字符串通过“-”连接，便于使用我们上述定义的P,R变量
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
        # 随机选择一个除了s5意外的状态s作为起点
        s = S[np.random.randint(4)]
        # 当前状态为终止状态或者时间步长太长时，一次采样结束
        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
            # 把(s,a,r,s_next)元组放入序列中
            episode.append((s,a,r,s_next))
            # s_next变成当前状态，开始接下来的循环
            s = s_next
        episodes.append(episode)
    return episodes

# 采样5次，每个序列最长不超过20步
episodes = sample(MDP,Pi_1,20,5)
print('第一条序列：\n',episodes[0])
print('第二条序列：\n',episodes[1])
print('第五条序列：\n',episodes[4])

---------------------------------------------------------------------
第一条序列：
 [('S1', '前往S2', 0, 'S2'), ('S2', '前往S1', -1, 'S1'), ('S1', '保持S1', -1, 'S1'), ('S1', '保持S1', -1, 'S1'), ('S1', '保持S1', -1, 'S1'), ('S1', '前往S2', 0, 'S2'), ('S2', '前往S3', -2, 'S3'), ('S3', '前往S4', -2, 'S4'), ('S4', '前往S5', 10, 'S5')]
第二条序列：
 [('S4', '概率前往', 1, 'S3'), ('S3', '前往S5', 0, 'S5')]
第五条序列：
 [('S2', '前往S1', -1, 'S1'), ('S1', '前往S2', 0, 'S2'), ('S2', '前往S3', -2, 'S3'), ('S3', '前往S5', 0, 'S5')]

计算所有状态的价值：

# 对所有采样序列计算所有状态的价值
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
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)

----------------------------------------------
使用蒙特卡洛法计算MDP状态价值为：
 {'S1': -1.186300691399171, 'S2': -1.6456669389656566, 'S3': 0.5300885668567163, 'S4': 5.916634388049788, 'S5': 0}