ホーム » 「ベイズ推定」タグがついた投稿

タグアーカイブ: ベイズ推定

バンディットアルゴリズムの続き

こんにちは。kzです。

今回は前回のバンディットアルゴリズムの続きです。UCBと簡単なトンプソンサンプリングの実装を行います。前回実装した\epsilon-GreedyとBoltzmann Softmaxとの最終的な比較も行います。

Keywords

この記事のキーワードを先にリストアップします。

  • \epsilon-Greedy
  • Boltzmann Softmax
  • Hoeffding’s inequality
  • UCB
  • Bayesian UCB
  • Thompson Sampling
  • scaled inverse chi-squared Distribution
先頭2つは前回やったので3つ目から話します。

Hoeffding’s inequality

Hoeffdingの不等式です。これはUCBの導出に使います。証明にはマルコフの不等式とチェビシェフの不等式を使います。まあそれは置いておいて、この不等式の意味を考えます。

確率変数Zの平均とその期待値の誤差t\in\mathbb{R}_+の確率をtで測っています。これがまさに後述するUCBの核の考え方です。こちら側でtそ操作することで、期待値のブレの許容範囲をある確率の範囲で操作することができます。

UCB

UCBはアームの期待値に加えてその上限(Upper Confidence Bound)、いわばアームの可能性(分散)、を考慮したアルゴリズムです。どうやって上限を加味するかといいますと、対象のアームを選択した回数の逆数で重み付けします。こうすることで多数引いたアームにはより小さい上限を、少数引いたアームにはより大きい上限を付与することができます。

これは直感的にもいいですよね。何度も引いたことのあるガチャガチャよりも数回しか引いたことのないガチャガチャの方が可能性を感じますよね。そういうことです。

ただし、
見ての通り、分母N_t(a)が小さいほど大きい値が加算されます。イメージは次の感じです。
お気づきかもしれませんが逆数と言いつつ、ルートを付けてます。気になりますよね。ということで証明します。先程のHoeffdingの不等式を用います。

Hoeffdingの不等式においてX_i \in [a,b][0,1]として次が得られます。(u\geq 0)。報酬はベルヌーイを仮定しました。
a \in \mathcal{A}を一つとり固定します。私たちの問題に書き換えるためNotationを定義します。

  • r_t(a) 報酬(確率変数)
  • \mathcal{Q}(a)  を真の行動価値
  • \mathcal{\hat{Q}}(a)  を経験上の行動価値
  • u をupper confidence bound, u=U_t(a)
代入しますと
ここで上限に対する閾値をこちらで設定します。
これをU_t(a)に関して解けば終了です
おや?よくみると2の位置が異なります。それは\epsilonをこちらで設定する必要があるからです。confidence levelと言われます。

UCBはこのconfidence levelをパラメータとして持ちます。もっとも一般的な値は次のようになります。
そしてこの時のアルゴリズムをUCB1いいます。代入して計算すると
先程の結果になりました。では実装してみます。
実はよくない結果です。前回のEGとBSと比べるとわかります。この記事の最後に全アルゴリズムの比較を行いますのでとりあえずは進めていきます。

Bayesian UCB

UCBにベイズを応用したアルゴリズムのことです。今までは報酬の分布にベルヌーイを仮定しましたが、他にもガウス分布の仮定があります。それぞれ次のように呼びます

  • Bernoulli – Thompson Sampling
  • Gaussian – Thompson Sampling
察しの通り事後分布からパラメータをサンプリングしてアームを選択します。

本当にハードなんですが、Gaussianの場合の詳しい計算は この方がやってくれています。

Thompson Sampling

定義は少し複雑ですが、実装の際はサンプリングを行うので思ったよりはシンプルです。
とあるaが他のどのa'よりも価値が高い、という事後確率(\pi)を計算してaを選ぶアルゴリズムです。R_tは時刻tまでのアームのhistoryです。

事後確率は積分を用いて次のように計算されます。
ここでR_iarm_iの報酬履歴です。左の積分が\mathcal{Q}(a) = \theta_iになる確率で右の積分がその他のアームで\mathcal{Q}(a') \leq \theta_iとなる確率です。積分が解析的ではないので、サンプリングします。困ったときはサンプリングです。結果アルゴリズムは次のように単純化できます。
  1. 各アームaの報酬の期待値に関する事後分布\pi( \theta  |R_a)から乱数\theta_aを生成
  2.  TS(a) = \arg \max_{a\in\mathcal{A}}  \theta_a \sim \pi( \theta |R_a)
何度も言いますが報酬の分布はベルヌーイを仮定するので事後分布であるベータ分布は非常にシンプルになります。
実装では一様分布を仮定するのでパラメータは共に1とします。

Best Arm RateをUCBと比較してみます。
全アルゴリズムでBARとCumulative Rewardsも比較してみます。
Thompson Samplingが圧勝でした。

scaled inverse chi-squared Distribution

もしGaussianでThompson Samplingを実装するならばこちらを\sigma^2の事前分布として導入します。(事前分布p(\mu, \sigma^2)を条件付き確率で分解する)scaled inverse chi-squared分布といいます。

    \[  p (\sigma^2) \sim  \mathcal{X}^2 (v_{0}, \hat{\sigma}_{0}^{2})  \]


Code

Gaussian Thompson Samplingすごく難しいですね、パラメータが多すぎてなにがなんだが、、、

次回はContextual BanditとIPSについてまとめたいと思います。でわ

References

MCMC入門 Gibbs Sampling

こんにちは。kzです。

本日はついにMCMCです。ベイズ推定とMAP推定がまだの方は次のリンクからどうぞ。

ギブスサンプリングをやっていくんですけど、その前にMCMCについて簡単に説明します。本当に簡単にです。

MCMC

マルコフ連鎖モンテカルロ法(Markov chain Monte Carlo methods、MCMC) の略であり、求めたい確率分布をマルコフ連鎖を作成することをもとに、確率分布のサンプリングを行うアルゴリズムの総称である。

ここでマルコフ連鎖とは
確率過程の一種であり、未来t+1の挙動が現在tの値だけで決定され、過去の挙動と無関係である(マルコフ性)。各時刻において起こる状態変化(遷移または推移)に関して、マルコフ連鎖は遷移確率が過去の状態によらず、現在の状態のみによる系列である。
難しそうに聞こえますが簡単です。明日の飯が今日の気分だけで決まる感じです。
ちなみにMCと確率遷移行列(マルコフ行列)と呼ばれる物はセットですが、ここでは割愛。でわでわ、気を取り直してMCMCです。

ベイズ推定について述べた際に、P(D)の部分が解析的に解けないことがよくあると言いました。たとば多変量になればなるほど積分がシンドくなりますよね。ということで確率分布を直接求めるのではなく、そこからサンプリングして擬似的に求めようというものなんです。

Gibbs Sampling

で、MCMCの一つであるギブスサンプリングを今回やります。これは確率分布を条件付き確率分布に無理やり分解して一変量ずつサンプリングしようという作戦。もちろん、マルコフ性を仮定してます。

たとえばp(x,y)という2変量の同時分布からサンプリングするよりもp(x|y)の方が簡単ですよね。次元が一つ減ってるようなものなので。で、このアルゴリズムはCoordinate descentに似てます。注目する変数以外を固定して計算していくので座標降下法のようにジグザグにサンプリングされます。で、もちろん初期値は適当に選ぶんですけど、バーンインというものがあります。これは初期値から一定時間までの間に得られたサンプルは捨てようという期間です。

ギブスサンプリングをいきなりベイズ推定で使うのは難しいのでまずは簡単な例でアルゴリズムを理解してからベイズ推定に応用しようと思います。ではいきましょう。

2次元ガウス分布を用いた例

二次元ガウス分布(等高線プロットか3Dプロット)からサンプリングしようと思います。つまり、p(x,y)から点をどんどん取ってくるイメージです。最終的に下を目指します。(赤い点がサンプリングした点)
まずはD次元ガウス分布の式を思い出します。
今回はD=2でやります。平均と分散は次のように設定します。
ただし、-1 < a < 1とします。逆行列と行列式も先に計算しておきます。
では同時分布p(x,y)をベクトル表記からスカラー表記に書き下します
で、目標の条件付き確率は
となり、先ほどの式を平方完成すると分母の周辺化が容易になるので
と、平方完成しました。あとはガッツリ積分していって条件付き確率を導出します。「なに分布にしたがうのかなー」と思いながら計算するといいです。

分母の積分をやっていくんですがガウス積分を使います。ガウス積分自体は極座標と重積分でいわゆるI^2計算で求まるのでここでは割愛します。
今、x_1以外は固定しているので定数とみなせます。
あとはガウス積分を用いて
周辺化も計算できたことなので条件付き確率を改めて計算すると
これは\mu= ax_2, \sigma =1の一次元ガウス分布ですね。つまり、正規分布にしたがうことがわかりました。よかったです。これでサンプリングできますね。次元も一つ落ちました。

ではコードを書いていく前に。アルゴリズムの流れの確認です。
  1. x_1, x_2 をそれぞれ初期化
  2. p(x_1 \mid x_2) よりx1を発生させる
  3. x_1 = x_1 で更新
  4. p(x_2 \mid x_1) よりx2を発生させる
  5. x_2 = x_2 で更新
  6. 2.〜5.を繰り返す
上の真の分布に対してギブスサンプリングで得られたものは
いいかんじですね。今回は平均もゼロですし、シンプルだったので簡単でしたが一般化するともう僕はシンドイです。。。

ベイス推定に応用

ではギブスサンプリングを用いてパラメータの事後分布をサンプリングで得たいと思います。対象となるパラメータは一次元ガウス分布の\mu, \sigmaです。 データ、尤度、事前分布を次のように設定します。
では\muをサンプリングするための条件付き確率の式を導出します。上の例のように正規分布になってくれるといいなーと思いながら導出します。
分母は関係ないので消し去りました。では\muに注目して計算します
よかったです。正規分布になりました。(4行目から5行目にかけては\muについての平方完成)

このノリで\sigmaをサンプリングするための条件付き確率の式を導出します。何度も言いますがどんな分布になるか想像しながら計算します。

先ほどと異なり、注目すべき変数は\sigmaです。となるとですね、
ここから例えば\sigmaでは平方完成できないですよね。となると正規分布には間違いなく従わないでしょう。んーーー、困った。

んーーーーーーーーーーー、と悩んでいても仕方ないので対数をとってみましょう。
ふむふむ、全くわからないです。僕は閃かなかったのですがガンマ関数がヒントになるらしいです。

とりあえず分母を無視して
ここで必殺技です。先ほど「んーーー」と止まった式を\tau=1 / \sigma^{2}と置き換えると(このタウは精度と呼ばれる)
となって、\alpha, \betaに注目すると以下がわかりました。
これすごいですよね、対数の式から分布を決定するなんてめっちゃ賢くないですか?まあこれをコードでかいてやってみると
うまくいったようです。よかったよかった。ただnumpyのガンマ関数のパラメータが逆数になっていることに注意してください。詳しくはコードにリンク貼ってます。
疲れましたあ。これ簡単な例だったからよかったですけど実際自分で全部モデリングしたらしんどそうですね、例えばパラメータが2から4になったら、、、、、もう眠れません。でわ

Reference