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

タグアーカイブ: ベイズ

事前分布に深入りする

2020年2月20日 / コメントする

こんにちは。kzです。

本日はprior distributionについて深入りしてみようと思います。事前分布って不思議ですよね。個人的には共役事前分布って誰が見つけたのかなー？と永遠思ってます。事前分布のスキルを高めるということはベイズの力を高めるということだと思います。なのでやっていきましょう。

Conjugate Prior

欲しいものはいつもposterior distributionです。

p(\theta)とp(\theta|X)が同じ分布になる時、そのpriorをconjugate priorといいます。共役事前分布です。しかし、実世界においてconjugate priorで十分か、といいますとそうではありません。なので現状はMCMCを使った次のようなアルゴリズムを使って数値計算するらしいです。

  • Stan
  • WinBUGS（Bayesian inference Using Gibbs Sampling）
  • OpenBUGS（Bayesian inference Using Gibbs Sampling）
  • JAGS（Just Another Gibbs Sampler）
  • PyMC（Hamiltonian Monte Carlo）
  • INLA（Integrated nested Laplace approximation）
とはいえ共役事前分布はもう皆さんご存知で、飽きていると思うので少し特殊なPriorを見ていこうと思います。

Improper Prior

事前分布が以下を満たす時、improper priorと言います。
たとえば、

  • Beta(0,0)
  • Uniform(-\infty, +\infty)
などです。priorがimproperであることが問題かどうかは人によってことなるようです。ベイジアンのLuca Rossiniという方はimproperでも問題ないと言っているようです。とはいえ、ベイズの定理の右辺が求まらないのであれば事後分布を得ることはできないのでもともこもありません。

Uninformative Prior

とはいえパラメータの事前分布なんてわからないのにどうすればいいんだ、、、という時に使われるのが無情報事前分布です。non-informative priorともいいます。情報がない、以外にも次の理由で使われたりします。

  • priorの事前情報がない
  • posteriorに変な影響を与えたくない
  • データのみの情報でposteriorを決定したい
一様分布はいい例じゃないでしょうか。そして他には以下のようなものがあります。
ここではJeffry’s Priorについて深入りしようと思います。

Jeffrey’s Prior

以下を満たす時、パラメータ\thetaのpriorが近似的にnon informative priorとなり（次で話します）、そのpriorをJeffreyの事前分布といいます。
見ての通り、フィッシャー情報量のルートになっています。また、このpriorはリパラメタライゼーションに対して不変です。これについては下で述べます。

とはいえ、これだけだとなにを意味してしているかよくわかりません。なので簡単な例、二項分布でJeffrey priorを求めてみます。
尤度は上のものを考えます。Jeffrey priorをp_J(\theta)と表すことにします。
よってフィッシャー情報量は
したがって
とこの尤度に対するJeffrey priorを求めることができました。しかもよく見るとこれはBeta(1/2,1/2)ですね。ついてでに分布も見ておきます。
ちなみにリパラメタライゼーションについては全単射な関数gを考えます。このとき\phi = g(\theta)という変数変換にたいしてp_J(\phi)は次の二通りの方法で得られます。
まずは変数変換後の尤度を用いて計算する方法。そして
ヤコビを使う方法です。たとえば
というリパラメタライゼーションを考える例は参考文献にあるので気になる方はどうぞ。

Jeffrey’s Priorの意味

さて、これを使う意味ですが個人的には２つあると思っています。
  • リパラメタライゼーションによる新しいモデル構築が容易な点
  • posteriorとのKLを最大ができる点
まず、一つ目は対象のパラメータに対して別のモデルを考える時にヤコビをかけるだけで済むのがいい点の一つだと思います。

そして二つ目、こちらが各だと思っています。Jeffreyを用いることでそのposteriorとのKLを最大化させることができます（データが十分多い時）。これによってposteriorはpriorには極力影響されない一方で、最大限にデータの影響のみを受けた分布になります。つまり、データをはっきりと表現していることになります。
実はこの話はBernardo reference priorというものと関係があるようです。そのパラメータ\thetaが一変数という条件がJeffrey Priorと一致するようです。

こういう意味でJeffrey priorはuninformativeと言われているようです。ついでにパラメータをいじった事後分布の結果を見てみましょう。
５種類の事前分布に対する事後分布の変化の図です。Jeffrey Priorは右から二つ目です。データの情報を最大限に説明する分布と言いましたが、まさにそうなんじゃないかなと思いました。まず注目すべきはy=0,n=1です。一回の試行で裏の出た時ですが、パラメータはグイっと左によっています。右の一様分布ではそうなりません。

次に注目したのがy=1,n=2です。\alpha=\beta=0.001の同じ条件と比べると事後分布が綺麗に半円を描いています。関数が滑らかであることはいろいろと都合がいいと思うのでこれもJeffreyの良さかなと思いました。

お気づきかと思いますが、\alpha=\beta=0.333と比べるとあまり差がないように思います。僕もそう思いもう少し調べてみました。
見ての通り、事後分布にそれほど大きな差はないと思いました。よって僕の中の結論としてはデータ数が非常に少ない場合においてJeffrey Priorは力を発揮するのでは？です。

僕的まとめ

この記事でreference priorとjeffrey priorの二つの新しい事前分布を紹介しましたが、共役に加えて最終的なベストなものは何なのか、というのが気になる点です。パラメータについて考慮するものがなければ共役事前分布またはreference priorがいいなと思いました。ここで、考慮するものとはたとえばデータのタイプであったり、物理的な何かであったり詳しいことは僕にはわかりませんが事前情報というよりは制約的な形で考慮するものなんだと思います。たとえば、画像処理における事前分布だと以下のようなものがあります。 そしてこのJeffrey Priorですが改めてまとめますと、これベルナドのrerefence priorの単変量バージョンに相当します。そしてこれはフィッシャー情報量から導出されます。これはつまり今あるデータから事前分布を構築するということを意味します。したがってよく使われる直感、inductive biasというべきでしょうか、の影響を極限まで希薄した事前分布を意味します。また、無情報事前分布と言われている理由はデータが十分に大きい場合、事後分布とのKLを最大化させるためです。この証明は参考文献にあります。最終的には場合分けで解が得られていました。

話が長くなりましたが今回個人的には学ぶことが多かったです。情報幾何に少し近づけた気がします。でわ

References

Variational Bayes 入門

2020年2月13日 / コメントする

こんにちは。kzです。

本日は変分ベイズやります。ちなみに呼び方はいくつかありまして、変文ベイズ、変分推論、Variational Approximation（変分近似）など。

変分とは？

そもそも変分ってなに！って感じですよね。実はこれ微分方程式とか汎関数とかが関わってます。

例えば点推定である最尤法はスコア関数（対数尤度）を微分して解を求めますよね。変分法はこれを拡張した概念です。

つまり、関数の関数の微分が変分です。
変分とは汎関数の微分です。詳しくは調べていただくとすぐわかりますが例えば最短曲線など求めるときはオイラーラグランジュ方程式と変分法を用いたりします。

変分ベイズとは？

事後分布p(\theta, X)を求める手法の一つです。MCMCと同じ感じ？と思った方もおられると思います。はい、MCMCも事後分布を求める手法です。

ではもういきなりですがVBとMCMCどっちがいいの？という疑問が生まれます。
ずばり、一番の違いは速度です。一般的にはVBの方が早く、計算の多いMCMCは遅いです。データXが多いときはVBの方がいいかもしれません。ただ上にあるようにVBは漸近的に分布を保証するわけではありません。

まずは変分下限の登場です。EMアルゴリズムの時と似てます。
zは潜在変数（パラメータ）です。H[z] = - \mathbb{E}_q[\log q(z)]は前回登場したShannon entropyですね。最後の式が変分下限です、Lと置くことにします。
せっかくなのでKLダイバージェンスとの関係も見てみます。
目的は事後分布p(z|X)を近似することです。なのでq(z)でKL最小化を目指します。式の最後で出てきたLは先程と同じ。つまり、KLの最小化は変分下限の最大化に等しいことがわかりました。

しかもよくみてください、KLって汎関数ですよね。

Mean-Field Approximation

変分ベイズでは平均場近似という重要なポイントがあります。これはq(z)に与える次の仮定です。

つまり、各潜在変数が独立であると仮定します。すると変分下限は次のように書き直せます。
計算を整理するために次を用います。対象の変数以外での期待値です。
これを用いて先程の式をさらに展開していきます。
iに注目していただければわかりますが、くくりあげています。第二項はiを除くjの項です。ラグランジュを用いて次の目的関数を得ます。
ここで変分法を使います。具体的にはq_j(z_i)に関してEuler-Lagrangeを用います。
ここで、Z_iは正規化定数と呼ばれるものです。以上からわかるようにq_iの最適化には他のすべてのq_jが必要になるのでイテラティブな更新になります。

一変数ガウスの例

次のような例を考えます。

変分ベイズを用いてパラメーラ\mu, \tauの分布を逐次的に求めます。まず以下の尤度をばらしておきます。
また平均場近似より次を仮定します。
では先程の計算を用いてq_\muを計算します。
\muに関して平方完成して以下を得ます。
同様の作業を次は\tauについて行います。
ガンマ分布に注意して
以上より
がわかりました。これより各期待値も簡単にもとまります。
したがって
を用いて更新すればいいことがわかります。
もう一つの例として混合ガウス分布をVBを使ったEMで求める例がありますが、非常に難しいので詳しくは以下をご覧ください。

応用

難しすぎて実装できませんでしたが変分ベイズを使った応用はたとえば自然言語処理におけるLDAがあります。とても難しい、、、、でわ

References

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

2020年2月7日 / コメントする

こんにちは。kzです。

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

Keywords

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

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

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_iはarm_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

ベイジアン入門　ベイズ推定

2019年12月12日 / コメントする

こんにちは。kzです。

前回、MAP推定までやりました。 しかしプロットもコードも書かなかったのでMAP推定をプロットするところからはじめます。復習がてらにね
この青いプロットは今回のテーマであるベイズ推定です。まぁそれは置いておいて、「尤度x事前分布」の最大化がMAP推定でした。プロットすると上図のように山が最も高くなる点になります。

MAP推定とベイズ推定の違いを簡単に

MAP推定は右辺の分子の最大化でした。一方でベイズ推定とは分母を切り捨てずちゃんと計算して左辺をちゃんと導出します。（なぜ切り捨ててよかったというとP(D)はパラメータ\thetaに依存しないからです。）

ということはベイズの方が難しいんですよ。（ちなみに最尤法、MAP推定は値を推定するので点推定とも呼ばれます。）

なぜベイズ推定したいのか？

これはあくまで僕なりの理解と見解なのでそれを踏まえて聞いて欲しいです。

前述通りMAP推定よりベイズ推定の方が難しいです。理由は分母の計算があるからです。この計算は積分の形になり、さらにいうと計算できなくて困っているくらいなのです（MCMCとかで対応）そこまでしてベイズ推定するメリットはなんなのか？なぜMAP推定ではダメなのか？

僕の理解だとまず、MAPがダメというわけではないです。ベイズ推定の結果がMAP推定と同値的なことになることはあります。
例えばベイズ推定後のパラメータの分布が上のようになったとしましょう。するとベイズ推定最大の利点の一つは「平均をとることもできる」という点だと思っています。説明すると、MAP推定のみの場合は得られる情報は真ん中の山の頂点を指す点、つまり\theta=10のみです。

一方で、ベイズ推定をすると事後分布の全体像がわかるため、例えば図においていい山が３つあるので平均をとって

    \[ \theta=\frac{6+10+17}{3}=11\]

のようにパラメータをより主観的に微調整できるのです。なぜ主観的という言葉を使ったかというとそもそもベイズの根幹に

直感・主観

を反映させるというアイデアがあるからです。

ベイズ推定

最も簡単なのでコイン投げを例に取りましょう。なぜ簡単かというとパラメータが一つで済むからです。というのは表が出る確率を\thetaとすると自動的に裏が出る確率が1-\thetaということです。

加えて計算も楽なのです。（僕は数学が苦手なので複雑な理論とかはできない。）ここではパラメータ\thetaの分布をベイズ推定で求めると同時に、その分布の変化を可視化することで理解を深めます。

パラメータの確率密度関数をf_n、コイン投げの結果x_nをH(head)、T(tale)で表し、確率変数をcと表すことにします。

例えば

    \[f_n(\theta) = f_n(\theta|c=H)\]

は表が出た時のベイズ推定で得た分布ということです。これを踏まえるとn+1回目の更新では次のようになります。

    \[f_{n+1}(\theta) = f_{n+1}(\theta|c=x_n) = \frac{p(c=x_n|\theta) f_n(\theta)}{p(c=x_n)} \]

ただし

    \[ p(c=x_n) = \int_0^1 p(c=x_n|\theta)f_n(\theta) d\theta \]

では事前分布を決めましょう。これは上式のf_1に相当します。事前分布とはパラメータに仮定する分布でした。今回のパラメータ\thetaはコインが表をだす確率です。僕なら事前分布に一様分布を仮定します。なぜなら表も裏も同じくらい出るだろうと思っているからです。

え？一様分布はいや？なら例えば次のような例を考えてみましょう。コイン持っていた人間が表がめっちゃ好きな人間なら事前分布は右に偏ったものになりますよね。つまり一様分布とはかけ離れ\theta=1付近で最大値をとるようないわばデルタ分布とかを仮定することになります。しかし、話がややこしくなるんです。なので一様分布で許してください。

つまり

    \[ f_1(\theta) = 1 \]


前回の記事を読んでくれた方は共役事前分布は？と思ったかもしれないがいったん忘れてください。なぜならとことん話が簡単になるからですよ。

残るは尤度、上式のp(c=x_n)のみです。尤度とはなんだった？それは試行の結果をパラメータを用いて表現したモデルでした

つまり、x_1=Hだったとするとその事象に対する尤度は

    \[ p(c=H|\theta) = \theta \]

となります、実にシンプルだろ？もし二項分布とかでやると計算がうわああってなります。（サイコロのようによりマルチなものをやろうとするとディリクレ分布とかが絡んできて、うわあああってなります）

ではベイズ推定してみよう。前述通りx_1=Hというデータが得られたとします。まずは分母つまり次の積分を計算しましょう。

    \[ p(c=x_1) = \int_0^1 p(c=H|\theta)f_n(\theta) d\theta \]

代入すると

    \[ p(c=x_1) = \int_0^1 \left( \theta \times 1 \right) d\theta = \frac{1}{2} \]

では分子も計算します

    \[f_{2}(\theta) = f_{2}(\theta|c=x_1) = \frac{p(c=H|\theta) f_1(\theta)}{1/2} \]

    \[f_{2}(\theta) = \frac{ \theta \times 1}{1/2} = 2\theta \]

もとまりました。少しまとめてみましょうか。
  1. （試行前）裏表同じくらいで出るよなー？
  2. （試行）Headが出た
  3. （試行後）まじ？まさか表出やすい？けどゆうてまだ一回目やしな
今この3にいるわけです。ではもう一度ベイズ推定してみましょう。
f_2はx_2というデータを用いたベイズ推定における事前分布となることに注意します。

二回目の試行も一回目と同じくx_2=Hだったとします。すると尤度は

    \[ p(c=H|\theta) \]

なのでこれを用いて先ほどと同様に積分と分布を計算します

    \[ p(c=x_2) = \int_0^1 \left( \theta \times f_2(\theta) \right) d\theta = \frac{2}{3} \]

    \[f_{3}(\theta) = f_{3}(\theta|c=x_2) = \frac{p(c=H|\theta) f_2(\theta)}{2/3} \]

    \[f_{3}(\theta) = \frac{ \theta \times 2\theta}{2/3} = 3\theta^2 \]

となりました。先ほど同様にまとめます。
  1. （試行前）裏表同じくらいで出るよなー？
  2. （試行）Headが出た
  3. （試行後）まじ？まさか表出やすい？けどゆうてまだ一回目やしな
  4. （試行）Headが出た
  5. （試行後）うせやろ？これ表しかでーへんパターンやろ
こんな気持ちですよね。この気持ちが分布（密度関数f_n）として現れているんです。というのはf_nに注目してみると
  1. f_1 = 1
  2. f_2 = 2\theta
  3. f_3 = 3\theta^2
とこれ単調増加具合が増えてますよね。可視化してみます。
まとめますとベイズ推定で得られる事後分布とはデータを知った後のパラメータに対する気持ちの分布なのです。今の例だと二回も連続で表が出たので表の出る確率が1に近いほど確率密度が高いのです。

もっとシュミレーションすると
ディリクレ分布でもやりたいんですけど、もう計算が酷くて、、、それよりmcmcしたいですよね。数値計算で積分求めるなんて素敵だと思いませんか？でわ