資源 | 跟著Sutton經典教材學強化學習中的蒙特卡羅方法（程式碼例項）

用動態規劃去解決強化學習的相關問題基本夠了，但還是有很多限制。比如，你知道現實世界問題的狀態轉移機率嗎？你能從任意狀態隨機開始嗎？你的MDP是有限的嗎？

好訊息是，蒙特卡羅方法能解決以上問題！蒙特卡羅是一種估計複雜的機率分佈的經典方法。本文部分內容取自Sutton的經典教材《強化學習》，並提供了額外的解釋和例子。

初探蒙特卡羅

蒙特卡羅模擬以摩納哥的著名賭場命名，因為機會和隨機結果是建模技術的核心，它們與輪盤賭，骰子和老虎機等遊戲非常相似。

相比於動態規劃，蒙特卡羅方法以一種全新的方式看待問題，它提出了這個問題：我需要從環境中拿走多少樣本去鑑別好的策略和壞的策略？

這一次，我們將再引入回報的概念，它是長期的預期收益：

有時，如果環節不收斂，那麼我們使用折扣因子：

我們將這些回報G_t和可能的A_t聯絡起來試圖推匯出：

根據大數定律，當N逼近∞時，我們能夠得到準確的期望。我們記i次模擬下標為i。

現在，如果這是一個馬爾科夫決策過程(MDP)（99%的強化學習問題都是），那麼我們知道它展現出了強馬爾科夫性質，也即：

有了這些，我們可以很容易推導得到這樣一個事實，即期望中的是完全無關的，從現在開始，我們將使G_s指示從某個狀態開始的回報（移動那個狀態到t=0）。

解決值函式的一種經典方式是對第一次s的發生的回報進行取樣，也叫首次訪問蒙特卡羅預測。那麼一個找到最優V的一個有效的演算法如下：

pi=init_pi()
returns=defaultdict(list)
for i in range(NUM_ITER):
    episode=generate_episode(pi)#(1)
    G=np.zeros(|S|)
    prev_reward=0
    for (state,reward) in reversed(episode):
        reward+=GAMMA*prev_reward
        #breaking up replaces s eventually,
        #so we get first-visit reward.
        G[s]=reward
        prev_reward=reward
    for state in STATES:
            returns[state].append(state)
        V={state:np.mean(ret) for state, ret in returns.items()}

另一種方法是每次訪問蒙特卡羅預測，也就是你在每個環節中每一次發生s的回報都進行取樣。在兩種情況下，估計均成平方收斂於期望。

在蒙特卡羅方法的背景下，策略迭代的核心問題是，正如我們之前說過的，如何確保探索和開採？

一種補救大狀態空間探索的方法是，明確我們從一個特定的狀態開始並採取特定的行動，對所有可能性採用輪循方式對它們的回報取樣。這假定我們可以從任何狀態出發，在每一環節的開始採取所有可能的行動，這在很多情況下不是一個合理的假設。然而對於像21點紙牌遊戲這樣的問題，這是完全合理的，這意味著我們可以很容易地解決我們的問題。

在以下程式碼中，我們只需要對我們之前的程式碼(1)做一個快速的補丁：

# Before(Start at some arbitrary s_0,a_0)
episode=generate_episode(pi)
# After(Start at some specifics s,a)
episode=generate_episode(pi,s,a)
# loop through s, a at every iteration.

線上策略ε-貪婪策略 

如果我們不能假設我們可以從任何狀態開始並採取任意行動那怎麼辦呢？好吧，那麼，只要我們不要太貪婪或者探索所有的狀態無窮次，我們仍然可以保證收斂，對嗎？

以上是線上策略方法的主要屬性之一，線上策略方法試圖去改善當前執行試驗的策略，與此同時，離線策略方法試圖改善一種不同於正在執行試驗的策略的策略。

說到這裡，我們要規則化“不要太貪婪”。一種簡答的方法是使用我們之前學過的k臂老虎機-ε-貪婪方法。回顧一下，我們以ε的機率從一個給定狀態下所有行動的均勻分佈中挑選，以1-ε的機率我們選argmax_tq(s,a)行動。

現在我們問：對於蒙特卡羅方法，這是否收斂到最優π_*？答案是它會收斂，但不是收斂到那個策略。

我們從q和一個ε-貪婪策略π(s)開始：

再一次，我們得到這個結論ε-貪婪策略，跟其他貪婪策略一樣，對於V_π執行單調的改進。如果我們回退到所有的時間步，那麼我們得到：

這就是我們想要的收斂性。

然而，我們需要去發現這個策略實際上收斂到什麼。顯然，即使最優策略是確定性的，我們的策略也被迫是隨機的，不能保證收斂到π*。然而，我們可以修訂我們的問題：

假定不是我們的策略保持以機率ε的隨機性一致選擇行動，而是環境不管我們的策略的規定隨機選擇一個行動，那麼，我們能夠確保一個最優解。證明的大綱在(1)中顯示，如果等式成立，那麼我們π=π，因此我們有Vπ=Vπ於環境，這個等式在隨機性下是最優的。

離線策略：重要性取樣

讓我們介紹一些新的術語！

• π是我們的目標策略。我們正努力最佳化它的預期回報。

• b是我們的行為策略。我們使b產π以後會用到的資料。

• π(a|s)>0⇒b(a|s)>0 ∀a∈A。這是收斂的概念。

離線策略方法通常有2個或者多個智慧體，其中一個智慧體產生另一個智慧體需要的資料，我們分別叫它們行為策略和目標策略。離線策略方法比線上策略方法更異想天開，就像神經網路之於線性模型更異想天開。離線策略方法往往更強大，其代價是產生更高的方差模型和較慢的收斂性。

現在，讓我們討論重要性取樣。

重要性取樣回答了這個問題：“給定E_π[G]，E_π[G]是什麼？”換句話說，你怎樣使用從b的取樣中獲得的資訊去決定π的期望結果。

一種你能想到的直觀方法就是：“如果b選擇a很多，π選a很多，那麼b的行為對於決π的行為是很重要的！”，相反：“如果b選擇a很多，π不曾選擇a，那麼b在a上的行為π在a上的行為沒有什麼重要性”，很有道理，對吧？

所以這差不多就知道重要性取樣的比率是什麼概念了。給定一個軌跡，在給定策略π的條件下，這條準確的軌跡發生的機率為：

π和b之間的比率是：

普通的重要性取樣

現在，有很多方法可以利用給我們一個E_π[G]的很好的估計。最基礎的方法是利用普通的重要性取樣。假定我們取樣了N個環節：

定義S的第一次到達時間為：

我們想要估計V_π(s),那麼我們可以使用經驗均值去透過首次訪問方法估計值函式：

當然，這可以很容易地推廣到每次訪問方法，但是我想用最簡單的形式來表達我的意思。這說明我們需要不同的方式來衡量每一環節的收益，因為對於π更容易發生的軌跡相比那些永遠不會發生的需要賦予更多的權重。

這種重要性取樣的方法是一種無偏估計量，但它存在極大的方差問題。假定第k個環節的重要性比率是1000，這個數很大，但是確實可能發生。這是不是意味著獎勵也要1000倍甚至更多呢？如果我們只有一個環節，我們的估計是準確的。從長遠來看，因為我們有一個乘法關係，導致比率可能會爆炸式增長或者消失。這對於估計來說，有點令人擔憂。

加權重要性取樣

為了減小方差，一種簡單直觀的方法就是減少估計的大小，透過除以重要比率的大小的總和（有點像柔性最大啟用函式）：

這叫做加權重要性取樣，它是一種有偏估計（偏差漸進趨於0），但是減小了方差。在此之前，我們能夠得到一個普通估計量的病態無界方差，但是這裡的每個元素的最大權值都是1，透過此限制了方差有界。Sutton建議，在實踐中，總是使用加權重要性取樣比較好。

增值實現

與許多其它取樣技術一樣，我們可以逐步實現它。假設我們使用上一節的加權重要性取樣方法，那麼我們可以得到一些如下形式的取樣演算法：

其中W_k可以是我們的權重。

我們想基於N_n來構造N_n+1，這是可行的。記C_n為，我們會以如下方式持續更新這個計算和：

C_n的更新規則非常明顯：

現在V_n是我們的值函式，但是一個非常相似的類比也可以應用到我們的行為Q_n。

當我們更新值函式的時候，我們也能更新我們的策略π，我們能夠用雖舊但是很好用argmax_tq(s,a)來更新π。

折扣意識重要性取樣

到目前為止，我們已經計算了回報，並取樣了回報得到了我們的估計。然而我們忽視了G的內部結構。它真的只是折扣獎勵的求和，我們未能將它納入到比率中ρ。折扣意識重要性取樣將γ建模為終止的機率。環節的機率在一些時間步t終止，因此必須是一個幾何分佈geo(γ)：

全部回報可以認為是對隨機變數R_t求期望：

可以構造一個任意的裂項求和如下：

以此類推，我們可以看到，令k從x處開始，那麼我們有γ^x

將上式代入G得到：

這將導致R_t項中的等效係數1，γ，γ²等。這就意味著，我們現在能夠分解G_t，並且在重要性取樣比率中使用折扣。

現在，回憶我們之前得到的：

如果我們擴充套件G，我們會有：

注意我們是怎樣在所有的回報中使用相同的比率的。一些回報，G_t.t+1，被整個軌跡的重要性比率相乘，這在模型假設：γ是終止機率下是不正確的。直觀上，我們想要給G_t.t+1P_t.t+1，這很容易：

啊，好多了！這樣，每個部分回報都有他們正確的比率，這極大解決了無界方差問題。

單個獎勵重要性取樣

另一種緩解p和它的方差問題的方式，我們可以將G分解為各個獎勵，然後做一些分析，讓我們研究一下P_t.T-1G_t.T：

對於每一項，我們有P_t.T-1γ^kR_t+k+1。擴充套件p，我們發現：

在沒有常數γ^k的情況下求期望：

記住E(AB)=E(A)E(B)當且僅當它們是獨立的。顯然根據馬爾科夫性質，任意π(A_i|S_i)和b(A_i|S_i)都是獨立於R_t+k+1，（如果i≥t+k+1），且

（b也一樣）。由此我們能夠將它們分離出來，從而得到：

這個式子看起來也許非常醜，但是我們發現：

所以我們可以完全忽略後半部分，從而得到：

這是什麼意思呢？我們完全可以用期望來表示最初的和：

這又一次將減少我們估計量的偏差。

Python中的線上策略模型

因為蒙特卡羅方法通常都是相似的結構。我在Python中建立了一個離散蒙特卡羅類，可以用來插入和執行。

程式碼下載：

https://github.com/OneRaynyDay/MonteCarloEngine

"""
General purpose Monte Carlo model for training on-policy methods.
"""
from copy import deepcopy
import numpy as np

class FiniteMCModel:
    def __init__(self, state_space, action_space, gamma=1.0, epsilon=0.1):
        """MCModel takes in state_space and action_space (finite) 
        Arguments
        ---------
        state_space: int OR list[observation], where observation is any hashable type from env's obs.
        action_space: int OR list[action], where action is any hashable type from env's actions.
        gamma: float, discounting factor.
        epsilon: float, epsilon-greedy parameter.
        If the parameter is an int, then we generate a list, and otherwise we generate a dictionary.
        >>> m = FiniteMCModel(2,3,epsilon=0)
        >>> m.Q
        [[0, 0, 0], [0, 0, 0]]
        >>> m.Q[0][1] = 1
        >>> m.Q
        [[0, 1, 0], [0, 0, 0]]
        >>> m.pi(1, 0)
        1
        >>> m.pi(1, 1)
        0
        >>> d = m.generate_returns([(0,0,0), (0,1,1), (1,0,1)])
        >>> assert(d == {(1, 0): 1, (0, 1): 2, (0, 0): 2})
        >>> m.choose_action(m.pi, 1)
        0
        """
        self.gamma = gamma
        self.epsilon = epsilon
        self.Q = None
        if isinstance(action_space, int):
            self.action_space = np.arange(action_space)
            actions = [0]*action_space
            # Action representation
            self._act_rep = "list"
        else:
            self.action_space = action_space
            actions = {k:0 for k in action_space}
            self._act_rep = "dict"
        if isinstance(state_space, int):
            self.state_space = np.arange(state_space)
            self.Q = [deepcopy(actions) for _ in range(state_space)]
        else:
            self.state_space = state_space
            self.Q = {k:deepcopy(actions) for k in state_space}
            
        # Frequency of state/action.
        self.Ql = deepcopy(self.Q)
    def pi(self, action, state):
        """pi(a,s,A,V) := pi(a|s)
        We take the argmax_a of Q(s,a).
        q[s] = [q(s,0), q(s,1), ...]
        """
        if self._act_rep == "list":
            if action == np.argmax(self.Q[state]):
                return 1
            return 0
        elif self._act_rep == "dict":
            if action == max(self.Q[state], key=self.Q[state].get):
                return 1
            return 0
    def b(self, action, state):
        """b(a,s,A) := b(a|s) 
        Sometimes you can only use a subset of the action space
        given the state.

        Randomly selects an action from a uniform distribution.
        """
        return self.epsilon/len(self.action_space) + (1-self.epsilon) * self.pi(action, state)
    def generate_returns(self, ep):
        """Backup on returns per time period in an epoch
        Arguments
        ---------
        ep: [(observation, action, reward)], an episode trajectory in chronological order.
        """
        G = {} # return on state
        C = 0 # cumulative reward
        for tpl in reversed(ep):
            observation, action, reward = tpl
            G[(observation, action)] = C = reward + self.gamma*C
        return G
    def choose_action(self, policy, state):
        """Uses specified policy to select an action randomly given the state.
        Arguments
        ---------
        policy: function, can be self.pi, or self.b, or another custom policy.
        state: observation of the environment.
        """
        probs = [policy(a, state) for a in self.action_space]
        return np.random.choice(self.action_space, p=probs)
    def update_Q(self, ep):
        """Performs a action-value update.
        Arguments
        ---------
        ep: [(observation, action, reward)], an episode trajectory in chronological order.
        """
        # Generate returns, return ratio
        G = self.generate_returns(ep)
        for s in G:
            state, action = s
            q = self.Q[state][action]
            self.Ql[state][action] += 1
            N = self.Ql[state][action]
            self.Q[state][action] = q * N/(N+1) + G[s]/(N+1)
    def score(self, env, policy, n_samples=1000):
        """Evaluates a specific policy with regards to the env.
        Arguments
        ---------
        env: an openai gym env, or anything that follows the api.
        policy: a function, could be self.pi, self.b, etc.
        """
        rewards = []
        for _ in range(n_samples):
            observation = env.reset()
            cum_rewards = 0
            while True:
                action = self.choose_action(policy, observation)
                observation, reward, done, _ = env.step(action)
                cum_rewards += reward
                if done:
                    rewards.append(cum_rewards)
                    break
        return np.mean(rewards)
if __name__ == "__main__":
    import doctest
    doctest.testmod()

如果你想在不同的庫中使用它，可以自己嘗試一下。

舉例：21點紙牌遊戲

在這個例子中，我們使用OpenAI的gym庫。在這裡，我們使用一個衰減的ε-貪婪策略去解決21點紙牌遊戲：

import gym
env = gym.make("Blackjack-v0")

# The typical imports
import gym
import numpy as np
import matplotlib.pyplot as plt
from mc import FiniteMCModel as MC

eps = 1000000
S = [(x, y, z) for x in range(4,22) for y in range(1,11) for z in [True,False]]
A = 2
m = MC(S, A, epsilon=1)
for i in range(1, eps+1):
    ep = []
    observation = env.reset()
    while True:
        # Choosing behavior policy
        action = m.choose_action(m.b, observation)

        # Run simulation
        next_observation, reward, done, _ = env.step(action)
        ep.append((observation, action, reward))
        observation = next_observation
        if done:
            break

    m.update_Q(ep)
    # Decaying epsilon, reach optimal policy
    m.epsilon = max((eps-i)/eps, 0.1)

print("Final expected returns : {}".format(m.score(env, m.pi, n_samples=10000)))

# plot a 3D wireframe like in the example mplot3d/wire3d_demo
X = np.arange(4, 21)
Y = np.arange(1, 10)
Z = np.array([np.array([m.Q[(x, y, False)][0] for x in X]) for y in Y])
X, Y = np.meshgrid(X, Y)

from mpl_toolkits.mplot3d.axes3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X, Y, Z, rstride=1, cstride=1)
ax.set_xlabel("Player's Hand")
ax.set_ylabel("Dealer's Hand")
ax.set_zlabel("Return")
plt.savefig("blackjackpolicy.png")
plt.show()

當沒有可用的A時，我們得到一個非常漂亮的圖形，（網路中Z為Flase）。

我也寫了一個模型的快速的離線策略版本，還沒有潤色過，因為我只想得到一個效能基準，這是結果：

Iterations: 100/1k/10k/100k/1million.
Tested on 10k samples for expected returns.
On-policy : greedy
-0.1636
-0.1063
-0.0648
-0.0458
-0.0312
On-policy : eps-greedy with eps=0.3
-0.2152
-0.1774
-0.1248
-0.1268
-0.1148
Off-policy weighted importance sampling:
-0.2393
-0.1347
-0.1176
-0.0813
-0.072

對程式碼的更改實際上非常小，因為正如我所說的，蒙特卡羅取樣是與環境無關的，我們修改了這部分程式碼（除去繪圖部分）：

# Before: Blackjack-v0
env = gym.make("CliffWalking-v0")
# Before: [(x, y, z) for x in range(4,22) for y in   range(1,11) for z in [True,False]]
S = 4*12
# Before: 2
A = 4

們執行gym庫得到了E_π[G]為-17.0。還不錯！懸崖行走問題是一幅一些區域為懸崖一些區域為平臺的地圖。如果你在平臺上行走，你的每一步將獲得-1的獎勵，如果你從懸崖上摔下來，將得到-100的獎勵，將你降落在懸崖上時，你會回到起點，不管地圖多大，每一個環節-17.0似乎是近乎最優的策略。

我們可以看到，蒙特卡羅方法對於計算任意行為和觀察空間具有詭異的機率分佈的任務的最優值函式和行為值是一種很好的技術。在未來，我們會考慮蒙特卡羅方法更好的變體，但是這也是強化學習基礎知識中的一塊偉大的基石。