探索アルゴリズム

reinforcement-learning

MCTS

Monte Carlo Tree Search

アルゴリズム

auto mcts(node_t const root_node)
{
  auto best_score = -std::numeric_limits<double>::infinity();
  auto best_seq   =  std::vector<node_t>();
  do {
    // 木を降る
    auto seq = std::vector<node_t>(1, root_node);
    traverse(seq);

    // 末端の状態からロールアウトする
    rollout(seq);

    // スコアを計算する
    auto score = cumulative_reward(seq);

    // ノードにスコアを伝播する
    for (auto const & node : seq) {
      update(node, score);
    }

    // 最適戦略を更新
    if (best_score < score) {
      best_score = std::move(score);
      best_seq   = std::move(seq  );
    }
  } while (has_remaining_time());

  return std::make_pair(best_score, std::move(best_seq));
}
auto traverse(std::vector<node_t> & seq)
{
  auto node = seq.back();
  while (node.expanded()) {
    // 多腕バンディット (mab; multi-armed bandit) の手法に基づいて
    // UCB1 値や Thompson sampling などで次のノードを選ぶ
    node = choose_next_mab(node);
    seq.push_back(node); // 展開する (末端も)
  }

  return seq;  
}
auto rollout(std::vector<node_t> & seq)
{
  while (seq.back().has_next()) {
    // 次の状態をランダムに選択する
    seq.push_back(choose_next_randomly(seq.back()));
  }
}

参考文献

  1. Cameron Browne, Edward Powley, Daniel Whitehouse, Simon Lucas, Peter I. Cowling, Stephen Tavener, Diego Perez, Spyridon Samothrakis, Simon Colton, et al.. 2012. A survey of Monte Carlo tree search methods. IEEE Transactions on Computational Intelligence and AI in Games, Vol. 4, No. 1, pp. 1-43.

Discounted UCB

割引 (Discounted) UCB

あらまし

多腕バンディットにおける UCB 値で過去の値を割り引く方法。最近プレイアウトした腕のほうが評価が正確だろうという仮定に基づいたヒューリスティック。

アルゴリズム

\textrm{DUCB}_{t,j} \equiv \frac{r_{t,j}^{D}}{n_{t,j}^{D}} + C \sqrt{\frac{\log s_{t}^{D}}{n_{t,j}^{D}}}

ここで

\begin{aligned}
n_{t+1,j}^{D} &\leftarrow \beta n_{t,j}^{D} + 1_{t+1,j} \\
r_{t+1,j}^{D} &\leftarrow \beta r_{t,j}^{D} + R_{t+1,j}
\end{aligned}

とする。$1_{t+1,j}$ は $t+1$ 回目のプレイアウトで $j$ が選択されたとき $1$ となり、それ以外のとき $0$ となる。また、$R_{t+1,j}$ は $t+1$ 回目のプレイアウトで $j$ が選択され、かつ勝利したとき $1$ となり、それ以外のとき $0$ となる。

ただし、$n_{0,j}^{D} = 0$ および $r_{0,j}^{D} = 0$ とする。

考察

割引操作は移動平均を計算しているのに近いようにみえる。移動平均するとバイアスがのるので補正したほうがよい気がするが、していない。

つまり

\begin{aligned}
E[n_{t,j}]
  &= E\left[\sum_{i=0}^{t} \beta^{t-i} 1_{t,j}\right] \\
  &= E\left[1_{t,j}\right] \sum_{i=1}^{t} \beta^{t-i} \\
  &= E\left[1_{t,j}\right] \frac{1 - \beta^{t}}{1 - \beta}
\end{aligned}

および

\begin{aligned}
E[r_{t,j}]
  &= E\left[\sum_{i=1}^{t} \beta^{t-i} r_{t,j}\right] \\
  &= E\left[r_{t,j}\right] \sum_{i=1}^{t} \beta^{t-i} \\
  &= E\left[r_{t,j}\right] \frac{1 - \beta^{t}}{1 - \beta}
\end{aligned}

となるので

\textrm{DUCB}_{t,j} \equiv \frac{r_{t,j}^{D}}{n_{t,j}^{D}} + C \sqrt{\frac{\log \frac{1 - \beta}{1 - \beta^{t}} s_{t}^{D}}{\frac{1 - \beta}{1 - \beta^{t}} n_{t,j}^{D}}}

とすべきではないか?

実装

final class DUCB {

	private static int COUNT = 100;
	private static double BETA = 0.9;

	private static double eval(double p) {
		return (Math.random() < p) ? 1.0 : 0.0;
	}

	private static int choose(double[] x, double[] n, int[] t) {
		double s = 0.0;
		for (int i = 0; i < n.length; i++) {
			double v = (1 - BETA) / (1 - Math.pow(BETA, t[i]));
			s += n[i] * v;
		}

		int    idxMax = -1;
		double ucbMax =  Double.NEGATIVE_INFINITY;
		for (int i = 0; i < n.length; i++) {
			double v = (1 - BETA) / (1 - Math.pow(BETA, t[i]));
			double ucb = x[i] / n[i] + Math.sqrt(2.0 * Math.log(s) / (n[i] * v));
			if (ucbMax < ucb) {
				ucbMax = ucb;
				idxMax = i;
			}
		}
		return idxMax;
	}

	public static void main(String[] args) {
		double[] p = {0.9, 0.1, 0.5, 0.3, 0.2};

		double[] x = new double[p.length];
		double[] n = new double[p.length];
		int   [] t = new int   [p.length];
		
		double regret = p[0] * COUNT;
		for (int i = 0; i < p.length; i++) {
			double e = eval(p[i]);
			regret -= e;
			x[i] = e;
			n[i] = 1.0;
			t[i] = 1;
		}
		for (int s = p.length; s < COUNT; s++) {
			int i = choose(x, n, t);
			double e = eval(p[i]);
			regret -= e;
			x[i] = BETA * x[i] + e;
			n[i] = BETA * n[i] + 1.0;
			t[i] ++;
		}
		
		System.out.printf("regret = %f", regret / COUNT);
	}
}

これだと

double ucb = x[i] / n[i] + Math.sqrt(2.0 * Math.log(n.length));

と置いたのと同じになる――ので元のままでよかったようだ。ただ試した限りでは、このほうが regret が小さくなる。

直近の結果が重視されるので失敗することが自動的な調整として働いているのかも知れない。

結局

double ucb = x[i] * (1 - BETA) / (1 - Math.pow(BETA, t[i]));

でいいということになる。チェビシェフの不等式から誤り率を $P$ として分散を考慮すると

double ucb = x[i] * (1 - BETA) / (1 - Math.pow(BETA, t[i]));
ucb += Math.sqrt(ucb * (1 - ucb) / P);

と補正できるから

final class DUCB {

	private static int COUNT = 100;
	private static double BETA = 0.9;
	private static double P = 0.1;

	private static double eval(double p) {
		return (Math.random() < p) ? 1.0 : 0.0;
	}

	private static int choose(double[] x, int[] t, boolean useVar) {
		int    idxMax = -1;
		double ucbMax =  Double.NEGATIVE_INFINITY;
		for (int i = 0; i < x.length; i++) {
			double ucb = x[i] * (1 - BETA) / (1 - Math.pow(BETA, t[i]));
			if (useVar) {
				ucb += Math.sqrt(ucb * (1 - ucb) / P);
			}
			if (ucbMax < ucb) {
				ucbMax = ucb;
				idxMax = i;
			}
		}
		return idxMax;
	}

	private static double ducb(double[] p, boolean useV) {
		double[] x = new double[p.length];
		int   [] t = new int   [p.length];
		
		double regret = p[0] * COUNT;
		for (int i = 0; i < p.length; i++) {
			double e = eval(p[i]);
			regret -= e;
			x[i] = e;
			t[i] = 1;
		}
		for (int s = p.length; s < COUNT; s++) {
			int i = choose(x, t, useV);
			double e = eval(p[i]);
			regret -= e;
			x[i] = BETA * x[i] + e;
			t[i] ++;
		}

		return regret / COUNT;
	} 

	public static void main(String[] args) {
		double[] p = {0.3, 0.1, 0.1, 0.1, 0.2};

		double r1 = 0.0;
		double r2 = 0.0;
		for (int n = 0; n < 1000; ++n) {
			r1 += ducb(p, false);
			r2 += ducb(p, true);
		}
		r1 /= 1000;
		r2 /= 1000;
		System.out.printf("%f - %f\n", r1, r2);
	}
}

となる。こちらだともう少し regret が小さくなる傾向になった。

参考文献

  1. Kocsis, L., Szepesvari, C.. 2006. Discounted UCB. Video Lecture. In: The Lectures of PASCAL Second Challenges Workshop. (link)

Accelerated UCT

加速 (accelerated) UCT

あらまし

割引 UCB の考え方を UCT に発展させた方法。

アルゴリズム

\textrm{AUCB}_{j} \equiv r_{j}^{A} + C \sqrt{\frac{\log s}{n_{j}}}

の式にもとづいてノードを選択する。右辺の第二項は通常の UCT (UCB) と同じ。

r_{j}^{A} \equiv
  \frac{\sum_{i \in \text{合法手}} v_{t,i} r_{i}^{A}}
       {\sum_{i \in \text{合法手}} v_{t,i}}

ここで $v_{0,j} = 0$ であり、

v_{t+1,j} \leftarrow \lambda v_{t,j} + 1_{t+1,j}

と更新していく。$1_{t+1,j}$ は $t+1$ 回目のプレイアウトで $j$ 番目の腕が選択されたとき $1$ で、それ以外のとき $0$ となる。

参考文献

  1. Junichi Hashimoto, Akihiro Kishimoto, Kazuki Yoshizoe, Kokolo Ikeda. 2012. Accelerated UCT and its application to two-player games. Advances in Computer Games, Lecture Notes in Computer Science Volume 7168, pp 1-12.

NMCS

Nested Monte Carlo Search

あらまし

ランダムシミュレーションや UCT (Upper Confidence bound applied to Trees) でよい性能が得られなかった Morpion Solitaire というパズルゲームで高い性能が報告され [1]、AMAF (All Move As First) との組み合わせで世界記録を達成した [2]。ぷよぷよにも応用されている [3]。

アルゴリズム

// NMCS を利用して思考
auto think(node_t const root_node, int const level)
{
  std::pair<move_t, score_t> best_move =
    std::make_pair(null_move, -std::numeric_limits<score_t>::infinity());

  do {
    auto const move = nmcs(root_node, level);
    if (best_move.first < move.first) {
      best_move = std::move(move);
    }
  } while (has_remaining_time());

  return best_move;
}
auto nmcs(node_t const node, int const level)
{
  std::pair<move_t, score_t> best_move =
    std::make_pair(null_move, -std::numeric_limits<score_t>::infinity());

  for (auto const& move : node.children()) {
    auto const move
      = (level <= 0)
      ? rollout(node[move])             // ロールアウトで評価する
      : nmcs   (node[move], level - 1); // 再帰的に評価する
    if (best_move.first < move.first) {
      best_move = std::move(move);
    }
  }

  return best_move;
}

参考文献

  1. Tristan Cazenave. 2009. Nested Monte-Carlo search. International Joint Conference on Artificial Intelligence 2009 (IJCAI 2009), pp. 456-461.
  2. 秋山晴彦, 小谷善行. 2010. Nested Monte-Carlo 探索の AMAF を用いた探索数調整による改良. 東京農工大学工学部情報工学科 東京農工大学大学院工学府 Technical Report, No. 7.
  3. 齋藤晃介, 三輪誠, 鶴岡慶雅, 近山隆. 2013. Nested Monte Carlo search のぷよぷよへの適用. ゲームプログラミングワークショップ2013論文集, pp. 134-137.

NRPA

Nested rollout policy adaptation

あらまし

ロールアウトのサンプリング分布を動的に調整するように NMCS を拡張した手法。

ロールアウトにおいて、状態ノード $\textrm{node}$ で行動 $a$ を起こす確率を $\exp(p_{\textrm{node},a})$ に比例させる。そしてその値を、ロールアウトごとに更新する。

アルゴリズム

重み付けを考えたロールアウトは下記のようになる。

auto rollout(node_t node)
{
  using std::begin;
  using std::end;

  assert(node);

  std::vector<int> seq;
  while (!node.empty()) {
    // 累積確率質量関数をつくる
    std::vector<double> w(node.size());
    w[0] = 0;
    for (auto const a : node.children()) {
      w[a+1] = w[a] + std::exp(p[node,a]);
    }

    // random() は [0,1) の乱数を返す関数
    auto const a = std::distance(begin(w)+1,
      std::lower_bound(begin(w), end(w), w.back() * random()));
    node = node[a];
    seq.push_back(a);
  }

   // return score and sequence
  return std::make_pair(evaluate(node), seq);
}
auto nrpa(int const level, node_t const node, policy_t & p)
{
  auto best_score = -std::numeric_limits<double>::infinity();
  auto best_seq   =  std::vector<int>();

  for (int i = 0; i < N; ++i) {
    auto const r
      = (level == 0)
      ? rollout(node)
      : nrpa(level - 1, node, p);
    if (best_score < r.first) {
      best_score = std::move(r.first );
      best_seq   = std::move(r.second);
    }

    // update policy by current best-seq
    adapt(node, p, best_seq);
  }

  return std::make_pair(best_score, best_seq);
}

これだと nexted になっていない気がする……

// 方策を更新
auto adapt(node_t node, policy_t & p, std::vector<int> const seq)
{
  auto const q = p;
  for (auto const a : seq) {
    p[node, a] += Alpha;

    auto rcpz = 0.0;
    for (auto const i : node) {
      rcpz = std::exp(q[node, i]);
    }
    rcpz = 1.0 / rcpz;

    for (auto const i : node) {
      p[node, i] -= Alpha * rcp * std::exp(q[node, i]);
    }

    node = node[a];
  }
}

ここで、NAlpha は適当な定数。

参考文献

  1. Chris Rosin. 2011. Nested rollout policy adaptation for Monte Carlo tree search. International Joint Conferences on Artificial Intelligence 2011 (IJCAI 2011), pp. 649-654.
  2. Tristan Cazenave, Fabien Teytaud. 2012. Application of the nested rollout policy adaptation algorithm to the traveling salesman problem with time windows. Learning and Intelligent Optimization, Lecture Notes in Computer Science 2012, pp 42-54.

NMCTS

Nested Monte Carlo tree search

アルゴリズム

auto nmcts(int const level = 1,
           std::vector<node_t> const &
           current_seq = std::vector(1, root_node))
{
  auto best_score = -std::numeric_limits<double>::infinity();
  auto best_seq   =  std::vector<node_t>();
  do {
    // 木を降る
    auto seq = current_seq;
    traverse(seq); // MCTS と同じ

    double score;
    if (level == 0) {
      // 末端の状態からロールアウトする
      rollout(seq); // MCTS と同じ

      // スコアを計算する
      score = cumulative_reward(seq);
    } else {
      std::tie(score, seq) = nmcts(level - 1, seq); // ※ここが MCTS と違う
    }

    // ノードにスコアを伝播する
    for (auto const & node : seq) {
      update(node, score);
    }

    // 最適戦略を更新
    if (best_score < score) {
      best_score = std::move(score);
      best_seq   = std::move(seq  );
    }
  } while (!time_limit(level));

  return std::make_pair(best_score, std::move(best_seq));
}

参考文献

  1. Hendrik Baier, and Mark H.M. Winands. 2012. Nested Monte-Carlo tree search for online planning in large MDPs. 20th European Conference on Artificial Intelligence (ECAI 2012), pp. 109-114.

BMCTS

Beam Monte-Carlo Tree Search

あらまし

ノードに対する訪問回数を数えておいて、同一深さのノードへの訪問回数の合計が一定数以上になったら、訪問回数が少ないノードを削除する。

アルゴリズム

auto bmcts(node_t const root_node)
{
  auto best_score = -std::numeric_limits<double>::infinity();
  auto best_seq   =  std::vector<node_t>();
  do {
    // 木を降る
    auto seq = std::vector<node_t>(1, root_node);
    traverse(seq); // ※このときノードを訪問した回数を数える

    // 末端の状態からロールアウトする
    rollout(seq);

    // スコアを計算する
    auto score = cumulative_reward(seq);

    // ノードにスコアを伝播する
    for (auto const & node : seq) {
      update(node, score); // ※このとき同じ深さの訪問回数の合計が
                           //   定数 B 以上なら上位 W ノード以外を木から削除する
    }

    // 最適戦略を更新
    if (best_score < score) {
      best_score = std::move(score);
      best_seq   = std::move(seq  );
    }
  } while (has_remaining_time());

  return std::make_pair(best_score, std::move(best_seq));
}

参考文献

  1. Hendrik Baier, Mark H. M. Winands. 2012. Beam Monte-Carlo tree search. IEEE Conference on Computational Intelligence and Games, pp. 227-233.

DNG-MCTS

Dirichlet Normal-Gamma based Monte-Carlo Tree Search

あらまし

マルコフ決定過程 (MDPs; Markov decision processes) で、状態の累積報酬を確率分布として表現して、Monte-Carlo tree search しながら、パラメータをベイジアン的に推定する。

方策 $\pi$ を持っているとする。ある状態 $s$ の報酬 $X_{s,\pi}$ を正規分布として、そこからある行動 $a$ を起こしたとき得られる報酬 $X_{s,a,\pi}$ を混合正規分布とする。正規分布のパラメータの事前分布を正規ガンマ分布として、混合正規分布のパラメータの事前分布をディリクレ分布とする。

ベイズ推定

状態 $s$ における報酬の分布を正規分布で表すので、そのパラメータの事前分布として正規ガンマ分布を使う。

正規ガンマ分布のパラメータ 報酬 $r$ を観測したときの更新式
$\alpha_{s}$ $\alpha_{s} \leftarrow \alpha_{s} + \frac{1}{2}$
$\beta_{s}$ $\beta_{s} \leftarrow \beta_{s} + \frac{\lambda_{s}}{\lambda_{s} + 1} \frac{(r - \mu_{s})^2}{2}$
$\mu_{s}$ $\mu_{s} \leftarrow \frac{\lambda_{s} + r}{\lambda_{s} + 1} \mu_{s}$
$\lambda_{s}$ $\lambda_{s} \leftarrow \lambda_{s} + 1$

状態 $s$ で行動 $a$ を起こしたときの報酬の分布を混合正規分布で表すので、その混合係数の事前分布としてディリクレ分布を使う。

ディリクレ分布のパラメータ 次の状態 $s’$ を観測したときの更新式
$\vec{\rho}_{s,a}$ $\rho_{s,a}^{(s’)} \leftarrow \rho_{s,a}^{(s’)} + 1$

アルゴリズム

auto online_planning(state_t s, tree_t t, int max_horizon)
{
  // すべての状態についてパラメータ mu[s], lambda[s], alpha[s], beta[s] を初期化する
  // すべての状態と行動についてパラメータベクトル rho[s,a] を初期化する
  init(mu, lambda, alpha, beta, rho)

  do {
    dng_mcts(s, t, max_horizon);
  } while (has_remaining_time());

  return thompson_sampling<false>(s, max_horizon);
}
auto dng_mcts(state_t s, tree_t t, int max_horizon)
{
  if (max_horizon == 0) {
    return 0;
  }

  if (!t.contain(node(s, max_horizon))) {
    t.emplace_add(s, max_horizon);
    // max_horizon ステップのロールアウトをプレイして報酬 r を推定する
    return r;
  }

  action_t a = thompson_sampling<true>(s, max_horizon);
  // 行動 a を起こして報酬 R(s,a) を観測して次の状態 s' に遷移する
  r = R(s,a) + gamma * dirichlet_normal_gamma_monte_carlo_tree_search(s_dash, t, max_horizon - 1);

  // 確率分布のパラメータの分布を更新
  alpha [s] += 0.5;                                               // 正規ガンマ分布のパラメータ
                                                                  //   ガンマ分布部分の形状母数
  beta  [s] += 0.5 * lambda[s] * (mu[s] - r)^2 / (lambda[s] + 1); //          尺度母数
  mu    [s]  =      (lambda[s] *  mu[s] + r)   / (lambda[s] + 1); //  正規分布部分の平均パラメータ
  lambda[s] ++;                                                   // サンプル数
  rho[s,a][s_dash]++;                                             // ディリクレ分布のパラメータ
}
template <bool sampling>
auto thompson_sampling(state_t s, int horizon)
{
  double best_q = -std::numeric_limits<double>::infinity();
  action_t best_a;
  for (action_t a : s) { // 状態 s でとれる行動 a すべてについて
    auto q = qvalue<sampling>(s, a, horizon);
    if (best_q < q) {
      best_q = q;
      best_a = a;
    }
  }
  return best_a;
}
template <bool sampling>
auto qvalue(state_t s, action_t a, int horizon) {
  double r = 0;
  double z = 0;
  for (state_t s_dash : s) { // 状態 s から行動 a で遷移しうるすべての状態 s' について
    if (sampling) { 
      w[s_dash] ~ Dir(rho[s,a]); // w[s'] を Dir(rho[s,a]) からサンプリング
    } else {
      w[s_dash] = rho[s,a][s_dash];
    }
    r += w[s_dash] * value<sampling>(s_dash, horizon-1);
    z += w[s_dash];
  }
  return R(s,a) + gamma * r / z;
}
template <bool sampling>
auto value(state_t s, int horizon)
{
  if ((horizon == 0) || s.terminal()) {
    return 0;
  }

  if (sampling) {
    // 正規ガンマ分布 NormGamma(mu[s], lambda[s], alpha[s], beta[s]) から (mu, tau) をサンプリング
    return mu;
  } else {
    return mu[s]
  }
}

参考文献

  1. Aijun Bai, Feng Wu, and Xiaoping Chen. 2013. Bayesian mixture modeling and inference based Thompson sampling in Monte-Carlo tree search. Proceedings of the Advances in Neural Information Processing Systems (NIPS-13), pp. 1646-1654.

D2NG-POMCP

Dirichlet-Dirichlet Normal-Gamma based Partially Observable Monte-Carlo Planning

あらまし

MDPs 用のアルゴリズムである DNG-MCTS を、部分観測マルコフ決定過程 (POMDPs; Partially observable Markov decision processes) 用に拡張した D2NG-POMCP を紹介する。

現在の状態 $s$ が信念の分布 $b$ によって確率的に与えられると考える。

ベイズモデリング

信念 $b$ において行動 $a$ を起こしたときの不変報酬を $X_{b,a}$ とする (信念も行動も離散で有限なので、二次元配列になっているイメージ。状態 $s$ は多項分布である信念 b に従う)。報酬も離散で有限 (ある特定の値しかとらない) $\Omega_{I} = {r_{1}, \ldots, r_{k}}$ として、$X_{b,a}$ が多項分布 $\textrm{Multinomial}(p_{1}, \ldots, p_{k})$ に従うとする。ここで $p_{i} = \sum_{s \in \Omega} b_{s} 1_{R_{s,a} = r_{i}}$ と仮定した。ただし $\sum_{i=1}^{k} p_{i} = 1$ である。ベイズ的に考えてパラメータ $p$ がディリクレ分布に従う $(p_{1}, \ldots, p_{k}) \sim \textrm{Dirichlet}(\Psi_{b,a})$ と考え、報酬 $r$ を観測したら $\psi_{b,a,r} \leftarrow \psi_{b,a,r} + 1$ と更新する。

信念 $b$ と状態 $s$ から政策 $\pi$ のもとで得られる累積の報酬を $X_{s,b,\pi}$ として、これが正規分布に従う $X_{s,b,\pi} \sim N(\mu_{s,b}, \frac{1}{\tau_{s,b}})$ とする。また、そのパラメータは正規ガンマ分布に従う $(\mu_{s,b}, \tau_{s,b}) \sim \textrm{NormalGamma}\left(\mu_{s,b,0}, \lambda_{s,b}, \alpha_{s,b}, \beta_{s,b}\right)$ として、報酬 $v$ を観測したら $\alpha_{s,b} \leftarrow \alpha_{s,b} + \frac{1}{2}$, $\beta_{s,b} \leftarrow \beta_{s,b} + \frac{1}{2} \frac{\lambda_{s,b} (v - \mu_{s,b,0})^2}{\lambda_{s,b} + 1}$, $\mu_{s,b,0} \leftarrow \frac{\lambda_{s,b} \mu_{s,b,0} + v}{\lambda_{s,b} + 1}$, $\lambda_{s,b} \leftarrow \lambda_{s,b} + 1$ と順に更新する。

信念 $b$ から政策 $\pi$ のもとで得られる累積報酬を $X_{b,\pi}$ として、これが混合正規分布に従う $X_{b,\pi} \sim \sum_{s \in \Omega_{s}} b_{s} N\left(\mu_{s,b}, \frac{1}{\tau_{s,b}}\right)$ とする。また、信念 $b$ から政策 $\pi$ のもとで、行動 $a$ を起こしたあとに得られる累積報酬 $X_{b,a,\pi}$ を $X_{b,a,\pi} = X_{b,a} + \gamma X_{b',\pi}$ とする。このとき $X_{b,a,\pi}$ の期待値は、$E\left[X_{b,a,\pi}\right] = E\left[X_{b,a}\right] + \gamma \sum_{o \in \Omega_{o}} E\left[X_{b',\pi}\right] \Omega_{o \mid b, a} 1_{b' = \xi_{b,a,o}}$ となる。ここで $\Omega(\cdot \mid b,a) \sim \textrm{Dirichlet}(\rho_{b,a})$ であり、$\rho_{b,a} = (\rho_{b,a,o_{1}}, \ldots, \rho_{b,a,o_{n}})$ として、$o$ を観測したら $\rho_{b,a,o} \leftarrow \rho_{b,a,o} + 1$ と更新する。

Thompson sampling

  1. パラメータをサンプリングする
    1. $\{w_{b,a,o}\} \sim \textrm{Dirichlet}(\rho_{b,a})$
    2. $\{w_{b,a,r}\} \sim \textrm{Dirichlet}(\psi_{b,a})$
    3. $\{\mu_{s',b'}\} \sim \textrm{NormalGamma}(\mu_{s',b',0}, \lambda_{s',b'}, \alpha_{s',b'}, \beta_{s',b'})$, where $b' = \xi_{b,a,o}$
  2. $\tilde{Q}$ をサンプリングして最大となる行動を選ぶ
\tilde{Q}(b,a) = \sum_{r \in \Omega_{I}} w_{b,a,r} r + \gamma \sum_{o \in \Omega_{o}} w_{b,a,o} 1_{b' = \psi_{b,a,o}} \sum_{s' \in \Omega_{s}} {b'}_{s'} \mu_{s',b'}

アルゴリズム

auto online_planning(history_t, tree_t & T)
{
  do {
    s = sample(P(h)); // 状態をサンプリング
    d2ng_pomcp(s, h, T, 0);
  } while (has_remaining_time());
  return thompson_sampling<false>(h, 0);
}
auto d2ng_pomcp(state_t const s, histroy_t h, tree_t & T, depth_t const d)
{
  if (d >= H) {
    return 0;
  }

  // node <h> is not in tree T
  if (!T.contains(h)) {
    // init parameters
    mu0[s,h], lambda[s,h], alpha[s,h], beta[s,h]
    rho[h,a], psi[h,a]

    T.add(h);

    return play_rollout(H - d);
  }

  auto const a = thompson_sampling<true>(h, d);

  // execute a by simulation, and
  // get next state s', observation o and reward i
  h_dash = hao;
  P(h_dash).add(s_dash);

  auto const r = i + gamma * d2ng_pomcp(s_dash, h_dash, T, d + 1);

  // update parameters
  alpha [s,h  ] += 0.5;
  beta  [s,h  ] += 0.5 * lambda[s,h] * (mu0[s,h] - r)^2 / (lambda[s,h] + 1);
  mu0   [s,h  ]  =      (lambda[s,h] *  mu0[s,h] + r)   / (lambda[s,h] + 1);
  lambda[s,h  ] += 1;
  rho   [h,a,o] += 1;
  pxi   [h,a,i] += 1;
}
template <boolean Sampling>
auto thompson_sampling(history_t const h, depth_t d)
{
  auto best_a = 0;
  auto best_q = -std::numeric_limits<double>::infinity();
  for (auto && a : Omega_A) {
    auto const q = qvalue<Sampling>(h, a, d);
    if (best_q < q) {
      best_q = q;
      best_a = a;
    }
  }
  return a;
}
template <bool Sampling>
auto qvalue(history_t const h, action_t const a, depth_t const d)
{
  using std::begin;
  using std::end;

  auto r = 0.0;

  // \gamma \sum_{o \in \Omega_{o}} w_{b,a,o}
  //   1_{b' = \psi_{b,a,o}}
  //   \sum_{s' \in \Omega_{s}} {b'}_{s'} \mu_{s',b'}
  {
    auto const rcpz
      = (Sampling)
      ? 1.0
      : 1.0 / std::accumulate(begin(rho[h,a]), end(rho[h,a]), 0.0);

    for (auto && o : Omega_o) {
      auto const w_o
        = Sampling
        ? sample(rho[h,a]   ) // w[o] をディリクレ分布からサンプリング
        :        rho[h,a][o];

      r += rcpz * w[o] * value<Sampling>(hao, d + 1, sampling);
    }
    r *= gamma;
  }

  // \sum_{r \in \Omega_{I}} w_{b,a,r} r
  {
    auto const rcpz
      = Sampling
      ? 1.0
      : 1.0 / std::accumulate(begin(psi[h,a]), end(psi[h,a]), 0.0);
    for (auto && r : Omega_I) {
      auto const w_r
        = Sampling
        ? sample(psi[h,a]   ) // w_r をディリクレ分布からサンプリング
        :        psi[h,a][r];

      r += rcpz * w_r * r;
    }
  }
  return r;
}
template <bool Sampling>
auto value(history_t const h, depth_t const d)
{
  if (d == H) {
    return 0;
  }

  auto v = 0.0;
  for (auto && s : P(h)) {
    auto const mu_s
      = Sampling
      ? // mu[s] を正規ガンマ分布からサンプリング
        sample(mu0[s,h], lambda[s,h], alpha[s,h], beta[s,h])
      : mu0[s,h];

    v += mu_s;
  }

  return v / std::size(P(h));
}

参考文献

  1. Aijun Bai, Feng Wu, Zongzhang Zhang, and Xiaoping Chen. 2014. Thompson sampling based Monte-Carlo planning in POMDPs. Proceedings of the 24th International Conference on Automated Planning and Scheduling (ICAPS 2014), pp. 21-26.