Introduction to Multi-Armed Bandits——04 Thompson Sampling[2]

参考资料

Russo D J, Van Roy B, Kazerouni A, et al. A tutorial on thompson sampling[J]. Foundations and Trends® in Machine Learning, 2018, 11(1): 1-96.
ts_tutorial

项目代码地址: https://github.com/yijunquan-afk/bandit-learning.git

一、General Thompson Sampling

Thompson Sampling 可以有效应用于 Bernoulli bandit 以外的一系列在线决策问题，我们现在考虑一个更普适的设置。

假设智能体(agent) 从集合 X\mathcal{X}X 选取一连串的动作(action) x1,x2,x3,⋯,x_1, x_2, x_3,\cdots,x1,x2,x3,⋯, 并应用于一个系统。
行动集可以是有限的，如Bernoulli bandit，也可以是无限的。
在应用动作 xtx_txt 之后，智能体观察到一个结果 yty_tyt，这是由系统根据条件概率 qθ(⋅∣xt)q_\theta(\cdot | x_t)qθ(⋅∣xt) 随机生成的。
智能体获得奖励 rt=r(yt)r_t=r(y_t)rt=r(yt)，其中 rrr 是一个已知的函数。智能体最初不知道 θ\thetaθ 的值，并使用先验分布 ppp 表示他的不确定性。

算法 4.1 和 4.2 以一种抽象的形式提出了贪心和TS方法，以适应这个更普适的设置。两者的区别在于它们生成模型参数 θ^\hat{\theta}θ^ 的方式。贪婪算法将 θ^\hat{\theta}θ^ 作为 θ\thetaθ 相对于分布 ppp 的期望值，而TS从 ppp 中抽取一个随机样本。然后，这两种算法都采取了使各自模型的预期奖励最大化的动作。

如果有一组有限的可能观测值yty_tyt，这个期望值由以下公式给出
Eqθ^[r(yt)∣xt=x]=∑oqθ^(o∣x)r(o)(4.1)\mathbb{E}_{q_{\hat{\theta}}}[r(y_t) | x_t = x] = \sum_o q_{\hat{\theta}}(o|x) r(o) \tag{4.1} Eqθ^[r(yt)∣xt=x]=o∑qθ^(o∣x)r(o)(4.1)

分布 ppp 根据已实现的观测值 y^t\hat{y}_ty^t 的条件进行更新。如果 θ\thetaθ 被限制在一个有限的集合中，这个条件分布可以被贝叶斯规则写成
Pp,q(θ=u∣xt,yt)=p(u)qu(yt∣xt)∑vp(v)qv(yt∣xt)(4.2)\mathbb{P}_{p, q}(\theta = u | x_t, y_t) = \frac{p(u) q_u(y_t | x_t)}{\sum_v p(v) q_v(y_t | x_t)}\tag{4.2} Pp,q(θ=u∣xt,yt)=∑vp(v)qv(yt∣xt)p(u)qu(yt∣xt)(4.2)

带有β先验的Bernoulli bandit是上述表述的一个特例。在这种特殊情况下，动作的集合是 X={1,…,K}\mathcal{X} = \{1,\ldots,K\}X={1,…,K}，只有奖励被观察到，所以 yt=rty_t = r_tyt=rt 。观察和奖励是由条件概率 qθ(1∣k)=θkq_\theta(1|k) = \theta_kqθ(1∣k)=θk and qθ(0∣k)=1−θkq_\theta(0|k) = 1-\theta_kqθ(0∣k)=1−θk 来模拟的。先验分布由向量 α\alphaα 和 β\betaβ 决定。其概率密度函数为：
p(θ)=∏k=1KΓ(α+β)Γ(αk)Γ(βk)θkαk−1(1−θk)βk−1p(\theta) = \prod_{k=1}^K \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha_k)\Gamma(\beta_k)} \theta_k^{\alpha_k-1} (1-\theta_k)^{\beta_k-1} p(θ)=k=1∏KΓ(αk)Γ(βk)Γ(α+β)θkαk−1(1−θk)βk−1

换句话说，在先验分布下，θ\thetaθ 的组成部分是独立的，并且是β分布的，参数为α\alphaα 和 β\betaβ。

对于这个问题，贪心算法（算法4.1）和TS算法（算法4.2}）在每个 ttt 的迭代中，对于 k∈{1,⋯,K}k \in \{1,\cdots,K\}k∈{1,⋯,K} 的后验参数以(αk,βk)(\alpha_k,\beta_k)(αk,βk) 开始。贪心算法将 θ^k\hat{\theta}_kθ^k 设定为期望，Ep[θk]=αk/(αk+βk)\mathbb{E}_p[\theta_k] = \alpha_k/(\alpha_k+\beta_k)Ep[θk]=αk/(αk+βk)，而TS从参数为 (αk,βk)(\alpha_k,\beta_k)(αk,βk) 的beta分布中随机抽取 θ^k\hat{\theta}_kθ^k 。然后，每个算法都会选择使 Eqθ^[r(yt)∣xt=x]=θ^x\mathbb{E}_{q_{\hat{\theta}}}[r(y_t) | x_t = x] = \hat{\theta}_xEqθ^[r(yt)∣xt=x]=θ^x 达到最大的动作 xxx。在应用选定的动作后，观察到奖励 rt=ytr_t = y_trt=yt ，并根据以下原则更新分布参数

(α,β)←(α+rt1xt,β+(1−rt)1xt)(\alpha, \beta) \leftarrow (\alpha + r_t {\bf 1}_{x_t}, \beta + (1-r_t) {\bf 1}_{x_t}) (α,β)←(α+rt1xt,β+(1−rt)1xt)

其中 1xt{\bf 1}_{x_t}1xt是一个分量 xtx_txt 等于 111，所有其他分量等于 000 的向量。

算法4.1和算法4.2还可以应用在更复杂的场景，让我们来看看最短路问题。

二、在线最短路问题

2.1 问题描述

一个人每天早上往返于家与单位之间。她想沿着平均遍历时间最少的路通勤，但她不确定不同路线的遍历时间。她该如何有效地学习，并在大量的遍历中尽量减少总遍历时间?

对于最短路问题，模型构建如下：

一个图 G=(V,E)G = (V, E)G=(V,E) , NNN 个点 V={1,…,N}V = \{1,\ldots, N\}V={1,…,N}，边 EEE 以及平均遍历时间(mean travel times) θ∈RN\theta \in \mathbb{R}^{N}θ∈RN。
点 111 是起点，点 NNN 是终点。
一个动作是一连串从起点到终点的不同边。
采取一个动作 xtx_txt 以后，对于每一个经过的边 e∈xte \in x_te∈xt，智能体观察到一个遍历时间 yt,ey_{t,e}yt,e，该遍历时间是从具有平均值 θe\theta_eθe 的分布中独立抽样得到的。
遍历过程中，智能体产生消耗 ∑e∈xtyt,e\sum_{e \in x_t} y_{t,e}∑e∈xtyt,e，于是奖励可以量化为 rt=−∑e∈xtyt,er_t = -\sum_{e \in x_t} y_{t,e}rt=−∑e∈xtyt,e。

2.2 Independent Travel Times

考虑先验分布：θe\theta_eθe 是独立的与log-Gaussian-distributed（参数为 μe\mu_eμe 和 σe2\sigma_e^2σe2）。也就是说，ln⁡(θe)∼N(μe,σe2)\ln(\theta_e) \sim N(\mu_e, \sigma_e^2)ln(θe)∼N(μe,σe2) 是高斯分布。因此，E[θe]=eμe+σe2/2\mathbb{E}[\theta_e] = e^{\mu_e + \sigma_e^2/2}E[θe]=eμe+σe2/2。此外，我们认为 yt,e∣θy_{t,e}|\thetayt,e∣θ 在各个边e∈Ee \in Ee∈E 上是独立的，并且是对数高斯分布，参数为ln⁡(θe)−σ~2/2\ln(\theta_e) - \tilde{\sigma}^2/2ln(θe)−σ~2/2 和 σ~2\tilde{\sigma}^2σ~2，因此 E[yt,e∣θe]=θe\mathbb{E}[y_{t,e}|\theta_e] = \theta_eE[yt,e∣θe]=θe。共轭特性适应于在观察 yt,ey_{t,e}yt,e 时更新θe\theta_eθe 分布:

(μe,σe2)←(1σe2μe+1σ~2(ln⁡(yt,e)+σ~22)1σe2+1σ~2,11σe2+1σ~2)(4.3)(\mu_e, \sigma_e^2) \leftarrow \left(\frac{\frac{1}{\sigma_e^2} \mu_e + \frac{1}{\tilde{\sigma}^2} \left(\ln(y_{t,e}) +\frac{\tilde{\sigma}^2}{2}\right)}{\frac{1}{\sigma_e^2} + \frac{1}{\tilde{\sigma}^2}}, \frac{1}{\frac{1}{\sigma_e^2} + \frac{1}{\tilde{\sigma}^2}}\right) \tag{4.3} (μe,σe2)←σe21+σ~21σe21μe+σ~21(ln(yt,e)+2σ~2),σe21+σ~211(4.3)

对于先验分布和后验分布属于同类别的分布，则先验与后验称为共轭分布，而先验分布被称为似然函数的共轭先验。

让我们用一个例子进行具体的说明：考虑到一个人每天从家到公司。假设可能的路径用图 G=(V,E)G = (V, E)G=(V,E) 表示，同时假设这个人知道每个边 e∈Ee \in Ee∈E 边的遍历距离 ded_ede，但是对平均遍历时间不确定。对他来说，构建一个期望等于遍历距离的先验是很自然的。对于对数高斯先验，可以通过设置 μe=ln⁡(de)−σe2/2\mu_e = \ln(d_e) -\sigma_e^2/2μe=ln(de)−σe2/2来实现。此时，E[θe]=de\mathbb{E}[\theta_e] = d_eE[θe]=de。参数 μe\mu_eμe 和 σe2\sigma_e^2σe2 也表达了不确定性的程度；一条边的平均遍历时间的先验方差是(eσe2−1)de2(e^{\sigma_e^2}-1) d_e^2(eσe2−1)de2。对数正态分布。

算法4.1和4.2可以应用在最短路问题上。在每个 ttt 的迭代中，对于 e∈Ee \in Ee∈E 的后验参数以(μe,σe)(\mu_e,\sigma_e)(μe,σe) 开始。贪心算法将 θ^e\hat{\theta}_eθ^e 设定为期望，Ep[θe]=eμe+σe2/2\mathbb{E}_p[\theta_e] = e^{\mu_e + \sigma_e^2/2}Ep[θe]=eμe+σe2/2，而TS从参数为 μe\mu_eμe 和 σe2\sigma_e^2σe2 的对数高斯分布中随机抽取 θ^e\hat{\theta}_eθ^e 。然后，每个算法都会选择使 Eqθ^[r(yt)∣xt=x]=−∑e∈xtθ^e\mathbb{E}_{q_{\hat{\theta}}}[r(y_t) | x_t = x] = -\sum_{e \in x_t} \hat{\theta}_eEqθ^[r(yt)∣xt=x]=−∑e∈xtθ^e 达到最大的动作 xxx。在应用选定的动作后，观察到输出 yty_tyt ，并根据公式4.3更新分布参数 (μe,σe2)(\mu_e, \sigma_e^2)(μe,σe2) 。

代码实现与分析

图使用binomial bridge的形式。有20层，所以从源头到目的地有184,756条路径。先验参数设置为 μe=−12\mu_e = -\frac{1}{2}μe=−21 和 σe2=1\sigma_e^2 = 1σe2=1 ，于是 E[θe]=1\mathbb{E}[\theta_e] = 1E[θe]=1 ，e∈Ee \in Ee∈E，以及条件分布参数为 σ~2=1\tilde{\sigma}^2 = 1σ~2=1。每个数据点代表了一万次独立模拟的平均值。

实现Dijkstra单源最短路算法：

from __future__ import generatorsclass priorityDictionary(dict):def __init__(self):'''Initialize priorityDictionary by creating binary heap ofpairs (value,key). Note that changing or removing a dict entrywill not remove the old pair from the heap until it is found bysmallest() or until the heap is rebuilt.'''self.__heap = []dict.__init__(self)def smallest(self):'''Find smallest item after removing deleted items from front ofheap.'''if len(self) == 0:raise IndexError("smallest of empty priorityDictionary")heap = self.__heapwhile heap[0][1] not in self or self[heap[0][1]] != heap[0][0]:lastItem = heap.pop()insertionPoint = 0while 1:smallChild = 2 * insertionPoint + 1if smallChild + 1 < len(heap) and \heap[smallChild] > heap[smallChild + 1]:smallChild += 1if smallChild >= len(heap) or lastItem <= heap[smallChild]:heap[insertionPoint] = lastItembreakheap[insertionPoint] = heap[smallChild]insertionPoint = smallChildreturn heap[0][1]def __iter__(self):'''Create destructive sorted iterator of priorityDictionary.'''def iterfn():while len(self) > 0:x = self.smallest()yield xdel self[x]return iterfn()def __setitem__(self, key, val):'''Change value stored in dictionary and add corresponding pairto heap. Rebuilds the heap if the number of deleted items getslarge, to avoid memory leakage.'''dict.__setitem__(self, key, val)heap = self.__heapif len(heap) > 2 * len(self):self.__heap = [(v, k) for k, v in self.items()]self.__heap.sort()# builtin sort probably faster than O(n)-time heapifyelse:newPair = (val, key)insertionPoint = len(heap)heap.append(None)while insertionPoint > 0 and \newPair < heap[(insertionPoint - 1) // 2]:heap[insertionPoint] = heap[(insertionPoint - 1) // 2]insertionPoint = (insertionPoint - 1) // 2heap[insertionPoint] = newPairdef setdefault(self, key, val):'''Reimplement setdefault to pass through our customized __setitem__.'''if key not in self:self[key] = valreturn self[key]def Dijkstra(G,start,end=None):D = {}  # dictionary of final distancesP = {}  # dictionary of predecessorsQ = priorityDictionary()   # est.dist. of non-final vert.Q[start] = 0for v in Q:D[v] = Q[v]if v == end: breakfor w in G[v]:vwLength = D[v] + G[v][w]if w in D:if vwLength < D[w]:raise ValueError("Dijkstra: found better path to already-final vertex")elif w not in Q or vwLength < Q[w]:Q[w] = vwLengthP[w] = vreturn (D,P)def shortestPath(G,start,end):D,P = Dijkstra(G,start,end)Path = []while 1:Path.append(end)if end == start: breakend = P[end]Path.reverse()return Path

再来看环境(environment)：

import numpy as np
from base.environment import Environment
from collections import defaultdict
from graph.dijkstra import Dijkstraclass IndependentBinomialBridge(Environment):"""Graph shortest path on a binomial bridge.The agent proceeds up/down for n_stages, but must end with equal ups/downs.e.g. (0, 0) - (1, 0) - (2, 0) for n_stages = 2\        /(1, 1)We label nodes (x, y) for x=0, 1, .., n_stages and y=0, .., y_limy_lim = x + 1 if x < n_stages / 2 and then decreases again appropriately."""def __init__(self, n_stages, mu0, sigma0, sigma_tilde=1.):"""graph[node1][node2] 表示node1 和 node2之间的边距Args:n_stages: 阶段数必须为偶数mu0: 独立分布的边的先验均值sigma0: 独立分布的边的先验标准差sigma_tilde: 标准差的观察噪声"""assert (n_stages % 2 == 0)self.n_stages = n_stagesself.mu0 = mu0self.sigma0 = sigma0self.sigma_tilde = sigma_tildeself.nodes = set()self.graph = defaultdict(dict)self.optimal_reward = None  # 当我们计算最短路径时填充self._create_graph()_ = self.get_shortest_path()def get_observation(self):"""这里的观察值是阶段数"""return self.n_stagesdef _get_width_bridge(self, x):"""在阶段x时计算bridge的宽度Args:x: 阶段数"""depth = x - 2 * np.maximum(x - self.n_stages / 2, 0) + 1return int(depth)def _create_graph(self):"""随机初始化图"""# 初始化结点for x in range(self.n_stages + 1):for y in range(self._get_width_bridge(x)):node = (x, y)self.nodes.add(node)# 添加边for x in range(self.n_stages + 1):for y in range(self._get_width_bridge(x)):node = (x, y)# 右上的结点right_up = (x + 1, y - 1)# 正右的结点right_equal = (x + 1, y)# 右下的结点right_down = (x + 1, y + 1)if right_down in self.nodes:distance = np.exp(# np.random.randn: 返回一个或一组服从标准正态分布的随机样本值。self.mu0 + self.sigma0 * np.random.randn())self.graph[node][right_down] = distanceif right_equal in self.nodes:distance = np.exp(self.mu0 + self.sigma0 * np.random.randn())self.graph[node][right_equal] = distanceif right_up in self.nodes and right_equal not in self.nodes:# 向上走distance = np.exp(self.mu0 + self.sigma0 * np.random.randn())self.graph[node][right_up] = distancedef overwrite_edge_length(self, edge_length):"""用确切的值覆盖原先的边长Args:edge_length: edge_length[start_node][end_node] = distance"""for start_node in edge_length:for end_node in edge_length[start_node]:self.graph[start_node][end_node] = edge_length[start_node][end_node]def get_shortest_path(self):"""找到最短路Returns:path: 遍历的结点列表"""start = (0, 0)end = (self.n_stages, 0)# 使用Dijkstra方法找到最短路final_distance, predecessor = Dijkstra(self.graph, start, end)path = []iter_node = endwhile True:path.append(iter_node)if iter_node == start:breakiter_node = predecessor[iter_node]path.reverse()# 更新最优奖励self.optimal_reward = -final_distance[end]return pathdef get_optimal_reward(self):return self.optimal_rewarddef get_expected_reward(self, path):"""给定一个路径，获取奖励Args:path: 结点列表Returns:expected_reward: -路长"""expected_distance = 0# zip()是Python的一个内建函数，它接受一系列可迭代的对象作为参数，# 将对象中对应的元素打包成一个个tuple（元组），# 然后返回由这些tuples组成的list（列表）。for start_node, end_node in zip(path, path[1:]):expected_distance += self.graph[start_node][end_node]return -expected_distancedef get_stochastic_reward(self, path):time_elapsed = defaultdict(dict)for start_node, end_node in zip(path, path[1:]):mean_time = self.graph[start_node][end_node]lognormal_mean = np.log(mean_time) - 0.5 * self.sigma_tilde**2stoch_time = np.exp(lognormal_mean + self.sigma_tilde * np.random.randn())time_elapsed[start_node][end_node] = stoch_timereturn time_elapsed

接下来设计智能体。

import copy
import randomclass IndependentBBEpsilonGreedy():"""Independent Binomial Bridge Epsilon Greedy"""def __init__(self, n_stages, mu0, sigma0, sigma_tilde, epsilon=0.0):"""An agent for graph bandits.Args:n_stages: binomial bridge的阶段数 (必须是偶数)mu0: 先验的平均值sigma0: 先验的标准差sigma_tilde: 观察值的噪声epsilon: 用于选择的参数"""assert (n_stages % 2 == 0)self.n_stages = n_stagesself.mu0 = mu0self.sigma0 = sigma0self.sigma_tilde = sigma_tildeself.epsilon = epsilon# 使用任意初始值设置内部环境self.internal_env = IndependentBinomialBridge(n_stages, mu0, sigma0)# 将边的后验保存为后验belief的元组(mean, std)self.posterior = copy.deepcopy(self.internal_env.graph)for start_node in self.posterior:for end_node in self.posterior[start_node]:self.posterior[start_node][end_node] = (mu0, sigma0)def get_posterior_mean(self):"""获得每条边的后验均值Returns:edge_length: edge_length[start_node][end_node] = distance"""edge_length = copy.deepcopy(self.posterior)for start_node in self.posterior:for end_node in self.posterior[start_node]:mean, std = self.posterior[start_node][end_node]edge_length[start_node][end_node] = np.exp(mean + 0.5 * std**2)return edge_lengthdef get_posterior_sample(self):"""获得每条边的后验抽样Return:edge_length: edge_length[start_node][end_node] = distance"""edge_length = copy.deepcopy(self.posterior)for start_node in self.posterior:for end_node in self.posterior[start_node]:mean, std = self.posterior[start_node][end_node]edge_length[start_node][end_node] = np.exp(mean + std * np.random.randn())return edge_lengthdef update_observation(self, observation, action, reward):"""为binomial bridge更新观察值.Args:observation: 阶段数action: 智能体选择的路(未使用)reward: reward[start_node][end_node] = stochastic_time"""assert observation == self.n_stagesfor start_node in reward:for end_node in reward[start_node]:y = reward[start_node][end_node]old_mean, old_std = self.posterior[start_node][end_node]# 转换精度，便于计算old_precision = 1. / (old_std**2)noise_precision = 1. / (self.sigma_tilde**2)new_precision = old_precision + noise_precisionnew_mean = (noise_precision * (np.log(y) + 0.5 /noise_precision) + old_precision * old_mean) / new_precisionnew_std = np.sqrt(1. / new_precision)# 更新后验值self.posterior[start_node][end_node] = (new_mean, new_std)def _pick_random_path(self):"""在bridge中完全随机地选择一条路径"""path = []start_node = (0, 0)while True:path += [start_node]if start_node == (self.n_stages, 0):breakstart_node = random.choice(list(self.posterior[start_node].keys()))return pathdef pick_action(self, observation):"""贪心地选择后验均值的最短路径"""if np.random.rand() < self.epsilon:path = self._pick_random_path()else:posterior_means = self.get_posterior_mean()self.internal_env.overwrite_edge_length(posterior_means)path = self.internal_env.get_shortest_path()return pathclass IndependentBBTS(IndependentBBEpsilonGreedy):"""Independent Binomial Bridge Thompson Sampling"""def pick_action(self, observation):"""从后验中抽样"""posterior_sample = self.get_posterior_sample()self.internal_env.overwrite_edge_length(posterior_sample)path = self.internal_env.get_shortest_path()return path

设置实验，不需要管action。

from base.experiment import BaseExperimentclass ExperimentNoAction(BaseExperiment):def run_step_maybe_log(self, t):# 观察环境，选择臂observation = self.environment.get_observation()action = self.agent.pick_action(observation)# 计算有用的值optimal_reward = self.environment.get_optimal_reward()expected_reward = self.environment.get_expected_reward(action)reward = self.environment.get_stochastic_reward(action)# 使用获得的奖励和选择的臂更新智能体self.agent.update_observation(observation, action, reward)# 记录需要的值instant_regret = optimal_reward - expected_rewardself.cum_optimal += optimal_rewardself.cum_regret += instant_regret# 环境进化（非平稳实验中才会用到）self.environment.advance(action, reward)self.data_dict = {'t': (t + 1),'instant_regret': instant_regret,'cum_regret': self.cum_regret, 'cum_optimal': self.cum_optimal,'unique_id': self.unique_id}self.results.append(self.data_dict)

跑一下书上的例图。这里的试验次数设置为200次，书中的是跑了5000次的结果。

import pandas as pd
import plotnine as ggdef generateTSAgent(n_steps, n_stages, mu0, sigma0, sigma_tilde, jobs):results = []for job_id in range(jobs):agent = IndependentBBTS(n_stages, mu0, sigma0, sigma_tilde)# 初始化环境，产生图env = IndependentBinomialBridge(n_stages, mu0, sigma0, sigma_tilde)experiment = ExperimentNoAction(agent, env, n_steps=n_steps, seed=job_id, unique_id=str(job_id))experiment.run_experiment()results.append(experiment.results)df_agent = (pd.concat(results)).assign(agent='TS')return df_agentdef generateEpsilonAgent(n_steps, n_stages, mu0, sigma0, sigma_tilde, jobs, epsilon=0):results = []for job_id in range(jobs):agent = IndependentBBEpsilonGreedy(n_stages, mu0, sigma0, sigma_tilde, epsilon)env = IndependentBinomialBridge(n_stages, mu0, sigma0, sigma_tilde)experiment = ExperimentNoAction(agent, env, n_steps=n_steps, seed=job_id, unique_id=str(job_id))experiment.run_experiment()results.append(experiment.results)df_agent = (pd.concat(results)).assign(agent='greedy-'+str(epsilon))return df_agentdef generateAgents():n_stages = 20n_steps = 500mu0 = -0.5sigma0 = 1sigma_tilde = 1N_JOBS = 200agents = []agents.append(generateEpsilonAgent(n_steps, n_stages, mu0, sigma0, sigma_tilde, N_JOBS))agents.append(generateEpsilonAgent(n_steps, n_stages, mu0,sigma0, sigma_tilde, N_JOBS, epsilon=0.01))agents.append(generateEpsilonAgent(n_steps, n_stages, mu0,sigma0, sigma_tilde, N_JOBS, epsilon=0.05))agents.append(generateEpsilonAgent(n_steps, n_stages, mu0,sigma0, sigma_tilde, N_JOBS, epsilon=0.1))agents.append(generateTSAgent(n_steps, n_stages,mu0, sigma0, sigma_tilde, N_JOBS))df_agents = pd.concat(agents)return df_agents   def plotCompare1():df_agents = generateAgents()plt_df = (df_agents.groupby(['t', 'agent']).agg({'instant_regret': np.mean}).reset_index())p = (gg.ggplot(plt_df)+ gg.aes('t', 'instant_regret', colour='agent')+ gg.geom_line(size=1.25, alpha=0.75)+ gg.xlab('time period (t)')+ gg.ylab('per-period regret')+ gg.scale_colour_brewer(name='agent', type='qual', palette='Set1'))print(p)def plotCompare2():df_agents = generateAgents()df_agents['cum_ratio'] = (df_agents.cum_optimal - df_agents.cum_regret) / df_agents.cum_optimalplt_df = (df_agents.groupby(['t', 'agent']).agg({'cum_ratio': np.mean}).reset_index())p = (gg.ggplot(plt_df)+ gg.aes('t', 'cum_ratio', colour='agent')+ gg.geom_line(size=1.25, alpha=0.75)+ gg.xlab('time period (t)')+ gg.ylab('Total distance / optimal')+ gg.scale_colour_brewer(name='agent', type='qual', palette='Set1')+ gg.aes(ymin=1)+ gg.geom_hline(yintercept=1, linetype='dashed', size=2, alpha=0.5))print(p)plotCompare1()
plotCompare2()

上图展示了应用贪心和TS算法处理最短路问题的结果。

悔值图表明，TS的性能很快就会收敛到最佳状态，而贪心算法却远非如此。我们通过改变 ϵ\epsilonϵ，来查看 epsilon−Greedyepsilon-Greedyepsilon−Greedy 算法的实验结果。对于每一次遍历。该算法会以 1−ϵ1-\epsilon1−ϵ 的概率选择一条由贪心算法产生的路径(exploit)，剩下则随机选择一条路径(explore)。虽然这种形式的探索是有帮助的，但图中显示，学习的速度远比TS慢。这是因为 epsilon−Greedyepsilon-Greedyepsilon−Greedy 在选择探索路径方面并不明智。TS会将探索努力引向信息丰富的路径，而不是完全随机的路径。

第二幅图表明，TS会以较快时间收敛到最优路径。对于 epsilon−Greedyepsilon-Greedyepsilon−Greedy 方法来说，情况就不是这样了。

2.3 Correlated Travel Times

在上例的基础上，让 θe\theta_eθe 是独立的和对数高斯分布的，参数为 μe\mu_eμe 和 σe2\sigma_e^2σe2，让观察值(observation)特征化为：
yt,e=ζt,eηtνt,ℓ(e)θey_{t,e} = \zeta_{t,e} \eta_t \nu_{t,\ell(e)} \theta_e yt,e=ζt,eηtνt,ℓ(e)θe
其中，每个 ζt,e\zeta_{t,e}ζt,e 代表与边 eee 相关的特殊因子，ηt\eta_tηt 代表所有边共同的因子，ℓ(e)\ell(e)ℓ(e)表示边eee是否位于binomial bridge的下半部分。而 νt,0\nu_{t,0}νt,0 和 νt,1\nu_{t,1}νt,1 分别代表对上半部和下半部的边共同有影响的因子。我们认为每个ζt,e\zeta_{t,e}ζt,e，ηt\eta_tηt， νt,0\nu_{t,0}νt,0 和 νt,1\nu_{t,1}νt,1都是独立的对数高斯分布，参数为 −σ~2/6-\tilde{\sigma}^2/6−σ~2/6 和 σ~2/3\tilde{\sigma}^2/3σ~2/3。参数ζt,e\zeta_{t,e}ζt,e, ηt\eta_tηt, νt,0\nu_{t,0}νt,0 和 νt,1\nu_{t,1}νt,1的分布是已知的，只有对应于每个边的参数 θe\theta_eθe 必须通过实验学习。鉴于这些参数，yt,e∣θy_{t,e} | \thetayt,e∣θ 的边缘分布与2.2完全相同，尽管联合分布 yt∣θy_t | \thetayt∣θ 不同。

共同因子诱发了binomial bridge中遍历时间之间的相关性。ηt\eta_tηt 代表了影响各地交通状况的随机事件，如当天的天气。而nut,0nu_{t,0}nut,0和nut,1nu_{t,1}nut,1则分别反映了只对一半的边的交通状况有影响的事件。各自反映了只对二项式桥的一半边缘的交通状况有影响的事件。尽管在先验条件下，边的平均遍历时间是独立的，但相关的观察结果会引起后验分布中的依赖性。

共轭性质（Conjugacy properties）再次促进了后验参数的高效更新，让 ϕ,zt∈RN\phi, z_t \in \mathbb{R}^Nϕ,zt∈RN 定义如下：
ϕe=ln⁡(θe)andzt,e={ln⁡(yt,e)if e∈xt0otherwise.\begin{equation*} \phi_e = \ln(\theta_e) \qquad \text{and} \qquad z_{t,e} =\begin{cases} \ln(y_{t,e}) \qquad & \text{if } e \in x_t \\ 0 \qquad & \text{otherwise.} \end{cases} \end{equation*} ϕe=ln(θe)andzt,e={ln(yt,e)0if e∈xtotherwise.

定义一个 ∣xt∣×∣xt∣|x_t| \times |x_t|∣xt∣×∣xt∣ 的协方差矩阵 Σ~\tilde{\Sigma}Σ~ ，其元素为:
Σ~e,e′={σ~2for e=e′2σ~2/3for e≠e′,ℓ(e)=ℓ(e′)σ~2/3otherwise,\begin{equation*} \tilde{\Sigma}_{e,e'} = \begin{cases} \tilde{\sigma}^2 \qquad & \text{for } e=e' \\ 2 \tilde{\sigma}^2/3 \qquad & \text{for } e \neq e', \ell(e) = \ell(e') \\ \tilde{\sigma}^2/3 \qquad & \text{otherwise,} \end{cases} \end{equation*} Σ~e,e′=⎩⎨⎧σ~22σ~2/3σ~2/3for e=e′for e=e′,ℓ(e)=ℓ(e′)otherwise,

其中e,e′∈xte,e' \in x_te,e′∈xt 。

N×NN \times NN×N的精度矩阵(concentration matrix):
C~e,e′={Σ~e,e′−1if e,e′∈xt0otherwise,\begin{equation*} \tilde{C}_{e,e'} = \begin{cases} \tilde{\Sigma}^{-1}_{e,e'} \qquad & \text{if } e, e' \in x_t\\ 0 \qquad & \text{otherwise,} \end{cases} \end{equation*} C~e,e′={Σ~e,e′−10if e,e′∈xtotherwise,
其中e,e′∈Ee,e' \in Ee,e′∈E 。

ϕ\phiϕ 的后验分布是高斯的，均值向量为 μμμ 和协方差矩阵为 ΣΣΣ，并根据以下公式更新：
(μ,Σ)←((Σ−1+C~)−1(Σ−1μ+C~zt),(Σ−1+C~)−1)\begin{align} (\mu, \Sigma) \leftarrow \left( \left(\Sigma^{-1} + \tilde{C} \right)^{-1} \left(\Sigma^{-1} \mu + \tilde{C} z_t\right), \left(\Sigma^{-1} + \tilde{C}\right)^{-1}\right) \tag{4.4} \end{align} (μ,Σ)←((Σ−1+C~)−1(Σ−1μ+C~zt),(Σ−1+C~)−1)(4.4)

TS算法也可以高效的计算方式应用：每个 ttt 的迭代从后验参数 μ∈ℜN\mu \in \Re^Nμ∈ℜN 和 Σ∈ℜN×N\Sigma\in \Re^{N\times N}Σ∈ℜN×N 开始。首先从均值为 μ\muμ 、协方差矩阵为 Σ\SigmaΣ 的高斯分布中抽出一个向量 θ^\hat{\theta}θ^ ，然后为每个e∈Ee \in Ee∈E 设定 hatthetae=phi^ehat{theta}_e = \hat{phi}_ehatthetae=phi^e，从而抽取样本 θ^\hat{\theta}θ^。选择一个动作 xxx 来最大化 Eqθ^[r(yt)∣xt=x]=−∑e∈xtθ^e\mathbb{E}_{q_{\hat{\theta}}}[r(y_t) | x_t = x] = -\sum_{e \in x_t} \hat{\theta}_eEqθ^[r(yt)∣xt=x]=−∑e∈xtθ^e ，使用Djikstra算法或其他算法。在应用选定的行动后，观察到结果yty_tyt，并根据公式（4.4）更新belief分布参数 (μ,Σ)(\mu, \Sigma)(μ,Σ) 。

代码实现与分析

单源最短路算法和上面的一致，先来看一下环境的不同：

class CorrelatedBinomialBridge(IndependentBinomialBridge):""" A Binomial Bridge with corrrelated elapsed time of each edge."""def is_in_lower_half(self, start_node, end_node):"""检查边缘start_node——>end_node是否位于桥的下半部分。"""start_depth = self._get_width_bridge(start_node[0])end_depth = self._get_width_bridge(end_node[0])if start_node[1] > start_depth / 2:return Trueelif start_node[1] < start_depth / 2:return Falseelse:return (start_depth<end_depth and end_node[1]==(start_node[1]+1)) \or (start_depth>end_depth and end_node[1]==start_node[1])def get_stochastic_reward(self, path):"""选择一条路，获得一个随机奖励.Args:path - list of list-like path of nodes from (0,0) to (n_stage, 0)Returns:time_elapsed - dict of dicts for elapsed time in each observed edge."""#shared factors:all_edges_factor = np.exp(-(self.sigma_tilde**2) / 6 +self.sigma_tilde * np.random.randn() / np.sqrt(3))upper_half_factor = np.exp(-(self.sigma_tilde**2) / 6 + self.sigma_tilde *np.random.randn() / np.sqrt(3))lower_half_factor = np.exp(-(self.sigma_tilde**2) / 6 + self.sigma_tilde *np.random.randn() / np.sqrt(3))time_elapsed = defaultdict(dict)for start_node, end_node in zip(path, path[1:]):mean_time = self.graph[start_node][end_node]idiosyncratic_factor = np.exp(-(self.sigma_tilde**2) / 6 +self.sigma_tilde * np.random.randn() / np.sqrt(3))if self.is_in_lower_half(start_node, end_node):stoch_time = lower_half_factor * all_edges_factor * idiosyncratic_factor * mean_timeelse:stoch_time = upper_half_factor * all_edges_factor * idiosyncratic_factor * mean_timetime_elapsed[start_node][end_node] = stoch_timereturn time_elapsed

接下来设计智能体：

import copy
import numpy as np
import numpy.linalg as nplafrom collections import defaultdict
from base.agent import Agent
from graph.env_graph_bandit import CorrelatedBinomialBridge_SMALL_NUMBER = 1e-10
###############################################################################
# Helper functions for correlated agentsdef _prepare_posterior_update_elements(observation, action, reward, num_edges, \edge2index, sigma_tilde, internal_env):"""生成用于相关BB问题后验更新的浓度矩阵Inputs:observation - 观察数 (= n_stages)action - 选择的动作，即一条路reward - 观察到的每条边的奖励, 字典的字典num_edges - 总的边数edge2index - 将每条边映射到一个唯一的下标sigma_tilde - 噪声internal_env - 内部环境Return:更新平均向量时使用的向量，更新协方差矩阵和平均向量时使用的浓度矩阵"""# 为每条边生成局部浓度矩阵和对数奖励log_rewards = np.zeros(num_edges)local_concentration = np.zeros((observation, observation))first_edge_counter = 0for start_node in reward:for end_node in reward[start_node]:log_rewards[edge2index[start_node][end_node]] = \np.log(reward[start_node][end_node])secod_edge_counter = 0for another_start_node in reward:for another_end_node in reward[another_start_node]:if first_edge_counter == secod_edge_counter:local_concentration[first_edge_counter,secod_edge_counter] \= sigma_tilde ** 2elif internal_env.is_in_lower_half(start_node, end_node) \== internal_env.is_in_lower_half(another_start_node, another_end_node):local_concentration[first_edge_counter, secod_edge_counter] \= 2 * (sigma_tilde ** 2) / 3else:local_concentration[first_edge_counter, secod_edge_counter] \= (sigma_tilde ** 2) / 3secod_edge_counter += 1first_edge_counter += 1# 求局部浓度矩阵的逆local_concentration_inv = npla.inv(local_concentration)# 生成浓度矩阵concentration = np.zeros((num_edges, num_edges))first_edge_counter = 0for start_node in reward:for end_node in reward[start_node]:secod_edge_counter = 0for another_start_node in reward:for another_end_node in reward[another_start_node]:concentration[edge2index[start_node][end_node] \,edge2index[another_start_node][another_end_node]] \= local_concentration_inv[first_edge_counter,secod_edge_counter]secod_edge_counter += 1first_edge_counter += 1return log_rewards, concentrationdef _update_posterior(posterior, log_rewards, concentration):"""更新后验参数Input:posterior - 后验参数的形式为(Mu, Sigma, Sigmainv)log_rewards - 对每条遍历边观察到的延迟的日志concentration - 根据新的观测计算出的浓度矩阵Return:updated parameters: Mu, Sigma, Sigmainv"""new_Sigma_inv = posterior[2] + concentrationnew_Sigma = npla.inv(new_Sigma_inv)new_Mu = new_Sigma.dot(posterior[2].dot(posterior[0]) +concentration.dot(log_rewards))return new_Mu, new_Sigma, new_Sigma_invdef _find_conditional_parameters(dim, S):"""给定一个维协方差矩阵S，返回一个包含用于计算每个组件的条件分布的元素的列表。"""Sig12Sig22inv = []cond_var = []for e in range(dim):S11 = copy.copy(S[e][e])S12 = S[e][:]S12 = np.delete(S12, e)S21 = S[e][:]S21 = np.delete(S21, e)S22 = S[:][:]S22 = np.delete(S22, e, 0)S22 = np.delete(S22, e, 1)S22inv = npla.inv(S22)S12S22inv = S12.dot(S22inv)Sig12Sig22inv.append(S12S22inv)cond_var.append(S11 - S12S22inv.dot(S21))return cond_var, Sig12Sig22inv##############################################################################class CorrelatedBBTS(Agent):"""Correlated Binomial Bridge Thompson Sampling"""def __init__(self, n_stages, mu0, sigma0, sigma_tilde, n_sweeps=10):"""An agent for graph bandits.Args:n_stages - number of stages of the binomial bridge (must be even)mu0 - prior meansigma0 - prior stddevsigma_tilde - noise on observationn_sweeps - number of sweeps, used only in Gibbs sampling"""assert (n_stages % 2 == 0)self.n_stages = n_stagesself.n_sweeps = n_sweeps# 使用任意初始值设置内部环境self.internal_env = CorrelatedBinomialBridge(n_stages, mu0, sigma0)# 保存一个映射(start_node,end_node)——>R以简化计算self.edge2index = defaultdict(dict)self.index2edge = defaultdict(dict)edge_counter = 0for start_node in self.internal_env.graph:for end_node in self.internal_env.graph[start_node]:self.edge2index[start_node][end_node] = edge_counterself.index2edge[edge_counter] = (start_node, end_node)edge_counter += 1# 保存所有边的数量self.num_edges = edge_counter# 先验参数self.Mu0 = np.array([mu0] * self.num_edges)self.Sigma0 = np.diag([sigma0**2] * self.num_edges)self.Sigma0inv = np.diag([(1 / sigma0)**2] * self.num_edges)self.sigma_tilde = sigma_tilde# p后验分布保存为包含平均向量、协方差矩阵及其逆的三重分布self.posterior = (self.Mu0, self.Sigma0, self.Sigma0inv)# boostrap版本中使用的附加参数self.concentration_history = []self.log_reward_history = []self.history_size = 0def get_posterior_mean(self):"""获得每条边的后验均值Return:edge_length - dict of dicts edge_length[start_node][end_node] = distance"""edge_length = copy.deepcopy(self.internal_env.graph)for start_node in edge_length:for end_node in edge_length[start_node]:edge_index = self.edge2index[start_node][end_node]mean = self.posterior[0][edge_index]var = self.posterior[0][edge_index, edge_index]edge_length[start_node][end_node] = np.exp(mean + 0.5 * var)return edge_lengthdef get_posterior_sample(self):"""获得每条边的后验抽样Return:edge_length - dict of dicts edge_length[start_node][end_node] = distance"""# flattened sampleflattened_sample = np.random.multivariate_normal(self.posterior[0],self.posterior[1])edge_length = copy.deepcopy(self.internal_env.graph)for start_node in edge_length:for end_node in edge_length[start_node]:edge_length[start_node][end_node] = \np.exp(flattened_sample[self.edge2index[start_node][end_node]])return edge_lengthdef update_observation(self, observation, action, reward):"""更新观察值Args:observation - number of stagesaction - path chosen by the agent (not used)reward - dict of dict reward[start_node][end_node] = stochastic_time"""assert (observation == self.n_stages)log_rewards, concentration = _prepare_posterior_update_elements(observation,\action, reward, self.num_edges, self.edge2index, self.sigma_tilde, \self.internal_env)# 更新联合分布的均值和方差矩阵new_Mu, new_Sigma, new_Sigma_inv = _update_posterior(self.posterior, \log_rewards, concentration)self.posterior = (new_Mu, new_Sigma, new_Sigma_inv)def pick_action(self, observation):"""Greedy shortest path wrt posterior sample."""posterior_sample = self.get_posterior_sample()self.internal_env.overwrite_edge_length(posterior_sample)path = self.internal_env.get_shortest_path()return path

实验环境和上面的一样，接下来跑一下：

import pandas as pd
import plotnine as ggdef generateIndependentBBTS(n_steps, n_stages, mu0, sigma0, sigma_tilde, jobs):results = []for job_id in range(jobs):agent = IndependentBBTS(n_stages, mu0, sigma0, sigma_tilde)# 初始化环境，产生图env = IndependentBinomialBridge(n_stages, mu0, sigma0, sigma_tilde)experiment = ExperimentNoAction(agent, env, n_steps=n_steps, seed=job_id, unique_id=str(job_id))experiment.run_experiment()results.append(experiment.results)df_agent = (pd.concat(results)).assign(agent='misspecified TS')return df_agentdef generateCorrelatedBBTS(n_steps, n_stages, mu0, sigma0, sigma_tilde, jobs):results = []for job_id in range(jobs):agent = CorrelatedBBTS(n_stages, mu0, sigma0, sigma_tilde)env = IndependentBinomialBridge(n_stages, mu0, sigma0, sigma_tilde)experiment = ExperimentNoAction(agent, env, n_steps=n_steps, seed=job_id, unique_id=str(job_id))experiment.run_experiment()results.append(experiment.results)df_agent = (pd.concat(results)).assign(agent='coherent TS')return df_agentdef generateAgents():n_stages = 20n_steps = 500mu0 = -0.5sigma0 = 1sigma_tilde = 1N_JOBS = 200agents = []agents.append(generateIndependentBBTS(n_steps, n_stages, mu0, sigma0, sigma_tilde, N_JOBS))agents.append(generateCorrelatedBBTS(n_steps, n_stages, mu0,sigma0, sigma_tilde, N_JOBS))df_agents = pd.concat(agents)return df_agentsdef plotCompare1():df_agents = generateAgents()plt_df = (df_agents.groupby(['t', 'agent']).agg({'instant_regret': np.mean}).reset_index())p = (gg.ggplot(plt_df)+ gg.aes('t', 'instant_regret', colour='agent')+ gg.geom_line(size=1.25, alpha=0.75)+ gg.xlab('time period (t)')+ gg.ylab('per-period regret')+ gg.scale_colour_brewer(name='agent', type='qual', palette='Set1'))print(p)def plotCompare2():df_agents = generateAgents()df_agents['cum_ratio'] = (df_agents.cum_optimal - df_agents.cum_regret) / df_agents.cum_optimalplt_df = (df_agents.groupby(['t', 'agent']).agg({'cum_ratio': np.mean}).reset_index())p = (gg.ggplot(plt_df)+ gg.aes('t', 'cum_ratio', colour='agent')+ gg.geom_line(size=1.25, alpha=0.75)+ gg.xlab('time period (t)')+ gg.ylab('Total distance / optimal')+ gg.scale_colour_brewer(name='agent', type='qual', palette='Set1')+ gg.aes(ymin=1)+ gg.geom_hline(yintercept=1, linetype='dashed', size=2, alpha=0.5))print(p)plotCompare1()
plotCompare2()

比较表明，由于考虑了边游历时间之间的相互依赖性，结果得到了实质性的改进。