Monte Carlo
本次代码实战使用“21点”游戏作为环境。首先导入必要的包并熟悉环境。
import gym
import numpy as np
from collections import defaultdict
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
env = gym.make('Blackjack-v0')
observation = env.reset()
print(env.action_space, env.observation_space, sep='\n')
Discrete(2)
Tuple(Discrete(32), Discrete(11), Discrete(2))
本次实现的是基于 ϵ − g r e e d y \epsilon-greedy ϵ−greedy 同轨策略的首次访问型MC算法,回忆其算法:
因此我们需要定义 $\epsilon-greedy \quad policy $ 。注意这里采用的是嵌套函数。
def make_epsilon_greedy_policy(Q_table, nA, epsilon):
def generate_policy(observation):
prob_A = np.ones(nA)*epsilon / nA
optimal_a = np.argmax(Q_table[observation])
prob_A[optimal_a] += (1.0 - epsilon)
return prob_A
return generate_policy
MC算法是逐幕进行的,所以我们要根据策略来生成一幕数据。
这里要注意:generate_policy
是一个函数即make_epsilon_greedy_policy
的返回值。generate_policy
的返回值是
π
\pi
π 。这里循环了1000次只是为了确保能获得完整的一幕。
def generate_one_episode(env, generate_policy):
trajectory = []
state = env.reset()
for i in range(1000):
Pi_table = generate_policy(state)
action = np.random.choice(np.arange(len(Pi_table)), p=Pi_table)
next_state, reward, done, _ = env.step(action)
trajectory.append((next_state, action, reward))
if done:
break
state = next_state
return trajectory
接下来是MC控制的主体部分,我们要循环足够多的次数使得价值函数收敛,每次循环都首先根据策略生成一幕样本序列,然后遍历每个“状态—价值”二元组,并用所有首次访问的回报的平均值作为估计.
这里要注意: Return
和Count
是字典,每个“状态—价值”二元组是一个key,该二元组每一幕的回报是它的value,随着越来越多的迭代,根据大数定律,它的平均值会收敛到它的期望值。并且在下一轮迭代生成另外一幕样本序列的时候,generate_policy
函数会根据Q_table更新!
def MC_control(env, iteration_times=500000, epsilon=0.1, discount_factor= 1.0):
Return, Count, Q_table = defaultdict(float), defaultdict(float), defaultdict(lambda: np.zeros(env.action_space.n))
policy = make_epsilon_greedy_policy(Q_table, env.action_space.n, epsilon)
for i in range(iteration_times):
trajectory = generate_one_episode(env, policy)
s_a_pairs = set([(x[0], x[1]) for x in trajectory])
for state, action in s_a_pairs:
s_a = (state, action)
first_visit_id = next(i for i, x in enumerate(trajectory) if x[0] == state and x[1] == action)
G = sum([x[2]*(discount_factor**i) for i, x in enumerate(trajectory[first_visit_id:])])
Return[s_a] += G
Count[s_a] += 1.
Q_table[state][action] = Return[s_a] / Count[s_a]
return policy, Q_table
接下来是将价值函数可视化。
def plot_value_function(Q_table):
x = np.arange(12, 21)
y = np.arange(1, 10)
X, Y = np.meshgrid(x, y)
Z_noace = np.apply_along_axis(lambda x: Q_table[(x[0], x[1], False)], 2, np.dstack([X, Y]))
Z_ace = np.apply_along_axis(lambda x: Q_table[(x[0], x[1], True)], 2, np.dstack([X, Y]))
def plot_surface(X, Y, Z, title):
fig = plt.figure(figsize=(20, 12))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1,cmap=matplotlib.cm.coolwarm, vmin=-1.0, vmax=1.0)
ax.set_xlabel('\nPlayer Sum', fontsize=25)
ax.set_ylabel('\nDealer Showing', fontsize=25)
ax.set_zlabel('\nValue', fontsize=25)
ax.set_title('\n\n' + title, fontsize=25)
ax.view_init(ax.elev, -120)
fig.colorbar(surf)
plt.savefig('../RL/%s.png'%title)
plt.show()
plot_surface(X, Y, Z_noace, '{}(not use Ace)'.format('optimal value function'))
plot_surface(X, Y, Z_ace, '{}(use Ace)'.format('optimal Value function'))
optimal_policy, Q_table = MC_control(env)
V = defaultdict(float)
for state, action_value in Q_table.items():
V[state] = np.max(action_value)
plot_value_function(V)
参考资料
《强化学习原理与Python实现》肖智清