モンテカルロ法2020年11月15日 | |
はじめにSutton and Barto 著『強化学習』(Reinforcement Learning An Introduction) を参考に、モンテカルロ法について見る。 環境
モンテカルロ法モンテカルロ法モンテカルロ法 (Mote Carlo method、MC 法) は、とりあえず走らせてみて評価するという方法である。数回の経験では環境の特定の情報しか得られないだろうが、大量に経験すれば環境全体のことがなんとなくわかるだろうという発想である。 ある方策のもとで試しに走ってみた状態・行動・報酬の系列をエピソード (episode) と呼ぶ。状態だけで見ると、エピソードは強化学習問題の状態空間の部分空間である。つまり、十分な学習にはそれだけの量と広さのエピソードが必要ということである。 モンテカルロ法ではマルコフ性を仮定していないので、非マルコフの問題にも使えるという特徴がある。 初回訪問 MC 法と逐一訪問 MC 法まず、方策がわかっているとして、これを評価する方法を考える。初期状態を適当に決めて、方策を使って走らせてみる。そのときの報酬系列から各ステップの収益が計算できるので、それを用いて状態価値を計算できる。1 度走っただけだと、たまたま通った状態しか評価されないので、満遍なく評価するために何度も走らせる。そうして平均をとったものを状態価値とすればよい。1 度の走査で同じ状態に複数回訪れる場合がある。最初に訪れた段階で評価する方法は初回訪問 MC 法と呼ばれる。すべての訪問で評価する場合は、逐一訪問 MC 法と呼ばれる。 モンテカルロ ES 法方策を更新するには、初期状態をランダムに決めて、初回訪問 MC 法で行動価値を評価し、グリーディ方策をとるということを繰り返せばよい。これは、開始点探索 (exploring starts) を含んだ方法ということで、モンテカルロ ES 法と呼ばれる。 方策オン型モンテカルロ法初期状態をランダムに決めることは、必要でない場合があるし、そもそも不可能な場合もある。といって、特定の開始点からグリーディに方策を決めると、局所的な解に落ち着いてしまうことがある。これを避けるため、グリーディ方策ではなく、確率 $\varepsilon$ だけ行動をランダムに決める $\varepsilon$ グリーディ方策を用いることが考えられる。このような方策は、すべての行動が選択される可能性があるということで $\varepsilon$ ソフト ($\varepsilon$-soft) 方策と呼ばれる。 具体的には、状態 $s$ に対する行動 $a$ を直接選ぶのではなくて、行動確率としての方策 $\pi(s,a)$ を以下のように更新する。 $$ a^* = \arg \max_a Q(s, a) $$ $$ \pi(s, a) = \begin{cases} 1 - \varepsilon + \varepsilon/n_a & (a = a^*) \\ \varepsilon/n_a & (a \ne a^*) \end{cases} $$ここで、$n_a$ は行動の個数である。グリーディでない行動が確率 $\varepsilon$ でそれぞれ $1/n_a$ の確率で選ばれるので、要するに $\varepsilon/n_a$ で選ばれるというわけである。 この方法は、エピソードを生成する方策と推定される方策が同じということで、方策オン型 (on-policy) と呼ばれる。一方、エピソードを生成する方策と推定される方策を分けて考える方法が、次の方策オフ型 (off-policy) である。 方策オフ型モンテカルロ法エピソードの生成に使う方策である挙動方策 (behavior policy) と、推定される方策である推定方策 (estimation policy) を分離することを考える。これにより、$\varepsilon$ ソフト方策を用いながら、そのランダム性が最適方策への改善を阻む状況を避けることができる。 そのためには、方策 $\pi'$ の経験から方策 $\pi$ を改善する必要がある。生成されたエピソード群について、状態 $s$ への $i$ 番目の初回訪問と、それに続いて起こる状態と行動の系列を考える。方策 $\pi$、$\pi'$ について、その系列が起こる確率をそれぞれ $p_i(s)$、$p'_i(s)$ とする。方策 $\pi'$ によって状態 $s$ から観測された収益を $R_i(s)$ とする。これらから $\pi$ についての状態価値を推定するには、$p_i(s)/p'_i(s)$ で重み付き平均を取る。 $$ V(s) = \frac{\sum_{i=1}^{n_s} \frac{p_i(s)}{p'_i(s)} R_i(s)}{\sum_{i=1}^{n_s} \frac{p_i(s)}{p'_i(s)}} $$ここで、$n_s$ はエピソード数である。$p_i(s)/p'_i(s)$ は次式で計算する。 $$ \frac{p_i(s_t)}{p'_i(s_t)} = \prod_{k=t}^{T_i(s)-1} \frac{\pi(s_k, a_k)}{\pi'(s_k, a_k)} $$ここで、$T_i(s)$ は $i$ 番目のエピソードの終了時刻である。 状態価値を反復計算の形に書き換える。$w_i = p_i(s)/p'_i(s)$ とすると、$n_s+1$ 番目の状態価値は次のように書ける。 $$ V_{n_s+1}(s) = V_{n_s}(s) + \frac{w_{n_s+1}}{W_{n_s+1}}\{R_{n_s+1} - V_{n_s}(s)\} $$ $$ W_{n_s+1} = W_{n_s} + w_{n_s+1} $$行動価値の場合は、同様に次のようにすればよい。 $$ Q_{n_s+1}(s, a) = Q_{n_s}(s, a) + \frac{w_{n_s+1}}{W_{n_s+1}}\{R_{n_s+1} - Q_{n_s}(s, a)\} $$この方法による方策の改善は、エピソードの末尾でグリーディな行動が続く系列で行われる。というのも、探索的行動が取られている場合、$\pi(s, a) = 0$ になる ($w = 0$ になってしまう) からである。この場合、グリーディな行動が取られているなら $\pi(s, a) = 1$ なので、$w$ の分子は結局 1 であり、$\pi'(s, a)$ だけ考えればよい。 この方法の問題点は、エピソード自体の系列が長くても、方策改善にはその末尾をちょろっと使うだけなので、たくさんのエピソードが必要で、計算に時間がかかるということである。 格子世界2 次元の格子世界 (gridworld) を動き回る問題を考える。動く方向は上下右左の 4 方向とする。 Python で問題を記述する。必要なモジュールを準備する。 import numpy as np from tqdm import tqdm 格子世界を実装する。
class Gridworld:
def __init__(self, nx, ny, goals):
self.nx = nx
self.ny = ny
self.goals = []
for g in goals:
self.goals.append(g[0] + g[1]*self.nx)
self.stop = self.goals.copy()
def get_states(self):
return np.arange(self.nx*self.ny)
def get_actions(self):
return np.arange(0, 4)
def step(self, s, a):
assert(a < 4)
if s in self.goals:
r = 0
else:
i = s % self.nx
j = s//self.nx
if a == 0: # up
j += 1
elif a == 1: # down
j -= 1
elif a == 2: # right
i += 1
elif a == 3: # left
i -=1
i = max(i, 0)
i = min(i, self.nx - 1)
j = max(j, 0)
j = min(j, self.ny - 1)
s = i + self.nx*j
r = -1
return s, r
def display_grid(self):
for j in range(self.ny-1, -1, -1):
print("{} ".format(j % 10), end="")
for i in range(self.nx):
s = i + j*self.nx
if s in self.goals:
print("G ", end="")
else:
print(". ", end="")
print()
print(" ", end="")
for i in range(self.nx):
print("{} ".format(i % 10), end="")
print()
def display_grid_value(self, v, d=6):
for j in range(self.ny-1, -1, -1):
print("{} ".format(j), end="")
for i in range(self.nx):
s = i + j*self.nx
print(("{:>%d.2f} " % d).format(v[s]), end="")
print()
print(" ", end="")
for i in range(self.nx):
print(("{:>%d} " % d).format(i), end="")
print()
def display_grid_action(self, pi):
for j in range(self.ny-1, -1, -1):
print("{} ".format(j % 10), end="")
for i in range(self.nx):
s = i + j*self.nx
if s in self.goals:
print("G ", end="")
else:
if len(pi.shape) == 1:
a = pi[s]
else:
a = np.argmax(pi[s, :])
if a == 0: # up
print("↑ ", end="")
elif a == 1: # down
print("↓ ", end="")
elif a == 2: # right
print("→ ", end="")
elif a == 3: # left
print("← ", end="")
print()
print(" ", end="")
for i in range(self.nx):
print("{} ".format(i % 10), end="")
print()
格子の大きさを nx x ny とし、格子の位置を (i, j) で表すと、状態変数は s = i + j*nx とする。左下を (0, 0) とする。報酬はひとつ動くたびに -1 とし、ゴールでは 0 とする。 ここでは 4 x 4 の格子とし、ゴールは (0, 3), (3, 0) の 2 つとする。 nx = 4 ny = 4 goals = [(0, ny-1), (nx-1, 0)] gw = Gridworld(nx, ny, goals) gw.display_grid() 3 G . . . 2 . . . . 1 . . . . 0 . . . G 0 1 2 3 初回訪問 MC 法行動はランダムとして、初回訪問 MC 法で状態価値を計算後、グリーディ方策を求める。
def gridworld_first_visit_mc(gw, episodes, steps):
states = gw.get_states()
actions = gw.get_actions()
n = len(states)
v = np.zeros(n)
count = np.zeros(n, dtype=int)
pi = np.zeros(n, dtype=int)
for episode in tqdm(range(episodes)):
ss = []
rs = []
s = np.random.choice(states)
for i in range(steps):
a = np.random.choice(actions)
ss.append(s)
s2, r = gw.step(s, a)
rs.append(r)
if s in gw.stop:
break
s = s2
ns = len(ss)
ret = np.zeros(ns)
ret[-1] = rs[-1]
for i in range(ns - 2, -1, -1):
ret[i] = ret[i+1] + rs[i]
uss = list(set(ss)) # unique
for s in uss:
for i in range(ns):
if ss[i] == s:
count[s] += 1
v[s] += (ret[i] - v[s])/count[s]
break
for s in states:
qs = np.zeros(len(actions))
for a in actions:
s2, r = gw.step(s, a)
qs[a] = r + v[s2]
pi[s] = actions[np.argmax(qs)]
return v, pi
実行。 np.random.seed(0) v_fvmc, pi_fvmc = gridworld_first_visit_mc(gw, 10000, 100) ここでは再現性確保のため、乱数シードを固定している。実行関数の数字の 1 つめはエピソード数、2 つめは計算ステップ数である。 状態価値の表示。 gw.display_grid_value(v_fvmc)
3 0.00 -13.76 -19.68 -21.60
2 -13.32 -17.58 -19.83 -19.54
1 -19.57 -19.80 -17.84 -14.02
0 -21.62 -19.73 -13.97 0.00
0 1 2 3
動的計画法で計算した状態価値は、以下の通りなので、まあまあの値である。
3 0.00 -14.00 -20.00 -22.00
2 -14.00 -18.00 -20.00 -20.00
1 -20.00 -20.00 -18.00 -14.00
0 -22.00 -20.00 -14.00 0.00
0 1 2 3
方策の表示。 gw.display_grid_action(pi_fvmc) 3 G ← ← ↓ 2 ↑ ← ← ↓ 1 ↑ ↑ ↓ ↓ 0 ↑ → → G 0 1 2 3 動的計画法で計算した方策は、以下の通りである。ただし、行動価値が同じ方向が複数ある場合でも、1 つだけを表示している。 3 G ← ← ↓ 2 ↑ ← ↓ ↓ 1 ↑ ↑ ↓ ↓ 0 ↑ → → G 0 1 2 3 モンテカルロ ES 法モンテカルロ ES 法で方策を求める。
def gridworld_mces(gw, episodes, steps):
states = gw.get_states()
actions = gw.get_actions()
n = len(states)
m = len(actions)
v = np.zeros(n)
q = np.zeros((n, m))
count = np.zeros((n, m), dtype=int)
pi = np.random.randint(m, size=n)
for episode in tqdm(range(episodes)):
sas = []
rs = []
s = np.random.choice(states)
for i in range(steps):
a = pi[s]
sas.append((s, a))
s2, r = gw.step(s, a)
rs.append(r)
if s in gw.stop:
break
s = s2
ns = len(sas)
ret = np.zeros(ns)
ret[-1] = rs[-1]
for i in range(ns - 2, -1, -1):
ret[i] = ret[i+1] + rs[i]
usas = list(set(sas)) # unique
for sa in usas:
for i in range(ns):
if sas[i] == sa:
s, a = sa
count[s, a] += 1
q[s, a] += (ret[i] - q[s, a])/count[s, a]
pi[s] = actions[np.argmax(q[s, :])]
break
for s in states:
v[s] = np.max(q[s, :])
return v, pi
実行。 np.random.seed(0) v_mces, pi_mces = gridworld_mces(gw, 10000, 5) 状態価値。 gw.display_grid_value(v_mces)
3 0.00 -1.00 -2.00 -3.00
2 -1.00 -2.00 -3.00 -2.00
1 -2.00 -3.00 -2.00 -1.00
0 -3.00 -4.00 -1.00 0.00
0 1 2 3
動的計画法による結果は、以下の通りである。
3 0.00 -1.00 -2.00 -3.00
2 -1.00 -2.00 -3.00 -2.00
1 -2.00 -3.00 -2.00 -1.00
0 -3.00 -2.00 -1.00 0.00
0 1 2 3
方策。 gw.display_grid_action(pi_mces) 3 G ← ← ↓ 2 ↑ ↑ ↑ ↓ 1 ↑ ← → ↓ 0 ↑ ↑ → G 0 1 2 3 動的計画法による結果は、以下の通りである。 3 G ← ← ↓ 2 ↑ ↑ ↑ ↓ 1 ↑ ↑ ↓ ↓ 0 ↑ → → G 0 1 2 3 位置 (1, 0) の状態価値の値が異なる結果になっている。これは、非効率な結果に収束してしまっていると考えられる。とはいえ、一応ゴールにはたどり着いてはいる。 方策オン型モンテカルロ法方策オン型モンテカルロ法で方策を求める。ここでは、初期状態は各状態を順番に選ぶものとした。
def gridworld_on_policy_mc(gw, episodes, steps, eps=0.):
states = gw.get_states()
actions = gw.get_actions()
n = len(states)
m = len(actions)
v = np.zeros(n)
q = np.zeros((n, m))
count = np.zeros((n, m), dtype=int)
pi = np.zeros((n, m)) + 1./m
si = 0
for episode in tqdm(range(episodes)):
sas = []
rs = []
s = states[si]
si += 1
if si == n:
si = 0
for i in range(steps):
a = np.random.choice(actions, p=pi[s, :])
sas.append((s, a))
s2, r = gw.step(s, a)
rs.append(r)
if s in gw.stop:
break
s = s2
ns = len(sas)
ret = np.zeros(ns)
ret[-1] = rs[-1]
for i in range(ns - 2, -1, -1):
ret[i] = ret[i+1] + rs[i]
usas = list(set(sas))
for sa in usas:
for i in range(ns):
if sas[i] == sa:
s, a = sa
count[s, a] += 1
q[s, a] += (ret[i] - q[s, a])/count[s, a]
a_max = actions[np.argmax(q[s, :])]
pi[s, :] = eps/m
pi[s, a_max] += 1. - eps
break
for s in states:
v[s] = np.max(q[s, :])
return v, pi
実行。 np.random.seed(0) v_onmc, pi_onmc = gridworld_on_policy_mc(gw, 10000, 5, 0.01) ここでは、$\varepsilon$ を 0.01 に設定している。 状態価値。 gw.display_grid_value(v_onmc)
3 0.00 -1.00 -2.02 -3.04
2 -1.00 -2.01 -3.01 -2.02
1 -2.01 -3.03 -2.01 -1.00
0 -3.02 -2.01 -1.00 0.00
0 1 2 3
方策。 gw.display_grid_action(pi_onmc) 3 G ← ← ← 2 ↑ ↑ ← ↓ 1 ↑ ← ↓ ↓ 0 → → → G 0 1 2 3 モンテカルロ ES 法と違い、問題ない結果となっている。$\varepsilon$ の探索が効いていると考えられる。 方策オフ型モンテカルロ法方策オフ型モンテカルロ法で方策を求める。ここでは、初期状態は各状態を順番に選ぶものとした。挙動方策として、推定方策からソフト方策を作っている。
def gridworld_off_policy_mc(gw, episodes, steps, eps=0.):
states = gw.get_states()
actions = gw.get_actions()
n = len(states)
m = len(actions)
v = np.zeros(n)
q = np.zeros((n, m))
c = np.zeros((n, m))
pi = np.random.randint(m, size=n)
si = 0
for episode in tqdm(range(episodes)):
pi_b = np.zeros((n, m))
for s in states:
pi_b[s, :] = eps/m
pi_b[s, pi[s]] += 1. - eps
sas = []
rs = []
s = states[si]
si += 1
if si == n:
si = 0
for i in range(steps):
a = np.random.choice(actions, p=pi_b[s, :])
sas.append((s, a))
s2, r = gw.step(s, a)
rs.append(r)
if s in gw.stop:
break
s = s2
ns = len(sas)
ret = 0.
w = 1.
for i in range(ns - 1, -1, -1):
ret += rs[i]
s, a = sas[i]
w /= pi_b[s, a]
c[s, a] += w
q[s, a] += (w/c[s, a])*(ret - q[s, a])
pi[s] = actions[np.argmax(q[s, :])]
if a != pi[s]:
break
for s in states:
v[s] = np.max(q[s, :])
return v, pi
エピソードの末尾を用いるため、方策改善のループではエピソードの末尾から逆に回している。方策が変更されたらループを抜ける。本来の式からすると、方策が変更されたときの w を使ってよいのか、という疑問がわく。ただ、その w を使ったから方策が変更されたとも言えるので...。w の更新をループ脱出のあとに置くと、計算がうまくいっていないように見えるケースがあった。 実行。 np.random.seed(0) v_offmc, pi_offmc = gridworld_off_policy_mc(gw, 20000, 5, 0.01) 状態価値。 gw.display_grid_value(v_offmc)
3 0.00 -1.00 -2.00 -2.75
2 -1.00 -2.00 -2.99 -2.00
1 -2.00 -2.85 -2.00 -1.00
0 -2.23 -2.00 -1.00 0.00
0 1 2 3
方策。 gw.display_grid_action(pi_offmc) 3 G ← ← ↓ 2 ↑ ↑ ↓ ↓ 1 ↑ ↑ ↓ ↓ 0 ↑ → → G 0 1 2 3 この問題に関しては、方策オン型よりも、方策オフ型のほうが収束が遅い (エピソード数が多くいる) 感じがする。 崖のある格子世界崖がある格子世界を考える。
class CliffGridworld(Gridworld):
def __init__(self, nx, ny, goals):
super().__init__(nx, ny, goals)
self.cliff = list(range(1, nx - 1))
self.stop += self.cliff
def step(self, s, a):
assert(a < 4)
if s in self.goals:
r = 0
elif s in self.cliff:
r = -100
else:
i = s % self.nx
j = s//self.nx
if a == 0: # up
j += 1
elif a == 1: # down
j -= 1
elif a == 2: # right
i += 1
elif a == 3: # left
i -=1
i = max(i, 0)
i = min(i, self.nx - 1)
j = max(j, 0)
j = min(j, self.ny - 1)
s = i + self.nx*j
r = -1
return s, r
def display_grid(self):
for j in range(self.ny-1, -1, -1):
print("{} ".format(j % 10), end="")
for i in range(self.nx):
s = i + j*self.nx
if s in self.goals:
print("G ", end="")
elif s in self.cliff:
print("C ", end="")
else:
print(". ", end="")
print()
print(" ", end="")
for i in range(self.nx):
print("{} ".format(i % 10), end="")
print()
def display_grid_action(self, pi):
for j in range(self.ny-1, -1, -1):
print("{} ".format(j % 10), end="")
for i in range(self.nx):
s = i + j*self.nx
if s in self.goals:
print("G ", end="")
elif s in self.cliff:
print("C ", end="")
else:
if len(pi.shape) == 1:
a = pi[s]
else:
a = np.argmax(pi[s, :])
if a == 0: # up
print("↑ ", end="")
elif a == 1: # down
print("↓ ", end="")
elif a == 2: # right
print("→ ", end="")
elif a == 3: # left
print("← ", end="")
print()
print(" ", end="")
for i in range(self.nx):
print("{} ".format(i % 10), end="")
print()
格子世界の左下をスタートとし、右下をゴールとして、その間を崖とする。移動中の報酬は -1 だが、崖の報酬は -100 で、崖に入ったらエピソードは終了するものとする。 nx = 12 ny = 4 goals = [(nx-1, 0)] gw = CliffGridworld(nx, ny, goals) gw.display_grid() 3 . . . . . . . . . . . . 2 . . . . . . . . . . . . 1 . . . . . . . . . . . . 0 . C C C C C C C C C C G 0 1 2 3 4 5 6 7 8 9 0 1 上の "C" が崖である。スタート地点は (0, 0) とするが、計算上は特にスタート地点を設けないものとする。 方策オン型モンテカルロ法方策オン型は、格子世界の関数をそのまま使えばよい。ここでは $\varepsilon = 0.1$ とする。まず、ステップ数を 100 として計算する。 np.random.seed(0) v_onmc, pi_onmc = gridworld_on_policy_mc(gw, 300000, 100, 0.1) 方策。 gw.display_grid_action(pi_onmc) 3 ↓ ← ↓ → → → → → → → ↓ ↓ 2 → → → → → → → → → → → ↓ 1 ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ → → ↓ 0 ↑ C C C C C C C C C C G 0 1 2 3 4 5 6 7 8 9 0 1 崖から少し離れたところを通っている。 ステップ数を 200 とした場合。 np.random.seed(0) v_onmc2, pi_onmc2 = gridworld_on_policy_mc(gw, 500000, 200, 0.1) 方策。 gw.display_grid_action(pi_onmc2) 3 ↓ → ↓ ↓ → → → → → ↓ ↓ ↓ 2 ↓ ↓ → → → → → → → → → ↓ 1 → → ↑ ↑ ↑ → ↑ ↑ ↑ → → ↓ 0 ↑ C C C C C C C C C C G 0 1 2 3 4 5 6 7 8 9 0 1 やはり、少し崖から離れ気味である。これは、探索的な行動によって崖に落ちやすいためと考えられる。 方策オフ型モンテカルロ法方策オフ型も、格子世界の関数をそのまま使えばよい。$\varepsilon = 0.1$ とする。まず、ステップ数を 100 として計算する。 np.random.seed(0) v_offmc, pi_offmc = gridworld_off_policy_mc(gw, 300000, 100, 0.1) 方策。 gw.display_grid_action(pi_offmc) 3 ↓ ↓ → ↓ ↓ ← ↓ ↓ ↓ → ↓ ↓ 2 ↓ ↓ ↓ → ↓ ← ↓ → ↓ ↓ ↓ ↓ 1 → → → → → → → → → → → ↓ 0 ↑ C C C C C C C C C C G 0 1 2 3 4 5 6 7 8 9 0 1 最短の経路を通っている。 ステップ数 200 の場合。 np.random.seed(0) v_offmc2, pi_offmc2 = gridworld_off_policy_mc(gw, 700000, 200, 0.1) 3 ↓ → ← ↓ ↓ → → ↓ → → → ↓ 2 → ↓ ↑ → ↓ ↓ ↓ → ↓ ↓ → ↓ 1 → → → → → → → → → → → ↓ 0 ↑ C C C C C C C C C C G 0 1 2 3 4 5 6 7 8 9 0 1 やはり最短の経路を通っている。ただし、かなりエピソード数が必要で、計算時間もかかる。方策オン型と違い、最短経路を見つけ出せているのは、挙動方策と推定方策を分けているおかげで、求める方策が探索的な行動の影響を受けずに済んでいるためと考えられる。 まとめモンテカルロ法は、動的計画法と違い、とにかく走ってみて学習する方法である。妥協案に落ち着かないようにしながら、十分に広く学ぶには、ソフト方策で探索的な行動を入れながら、たくさんの試行を行う必要がある。学習の方法として、経験から直接学ぶ方策オン型と、経験から間接的に学ぶ方策オフ型がある。崖のある格子世界において、方策オン型だと安全を見た結果になったが、方策オフ型だと最適な結果が導かれるのを見た。例えるなら、方策オフ型だと痛い目を見る人と学ぶ人が違うので、学習に際して怖れる必要がないというわけである。ただ、方策オフ型は、収束が遅くエピソード数が多く必要なようである。 | |
| PENGUINITIS | |