トンプソン抽出の例で見る共役事前分布

※この記事は、Speee Advent Calendar 4日目の記事です。 昨日の記事はこちら

tech.speee.jp

初めまして。アドプラットフォーム事業部の本田 @mov_vc です。SpeeeでUZOUの開発をしています。

以前に下記の記事で紹介されていますが、UZOUで利用されているバンディットアルゴリズムはトンプソン抽出 (Thompson Sampling) で実装されています。

tech.speee.jp

大まかに言うと、トンプソン抽出によるバンディットアルゴリズムでは

1. それぞれのアームの期待値の確率分布に従って乱数を生成する
2. 生成された乱数が最大のアームを1回引いて期待値の確率分布を修正する

という手順を繰り返します。この記事では、これらがどのようにして実現可能なのかという理論面に焦点を当てて観察していきます。

この記事で話さないこと

• バンディットアルゴリズム自体の説明（先の記事参照）
• リグレット最小化の根拠（いつか話したい😢）

アルゴリズムを実現するために必要なこと

トンプソン抽出のアルゴリズムは、

1. それぞれのアームの期待値の確率分布に従って乱数を生成する
2. 生成された乱数が最大のアームを1回引いて期待値の確率分布を修正する

という2つの手順を繰り返すものです。つまり、アルゴリズムを実現するために必要なのは、

• それぞれのアームの期待値の確率分布
• 期待値の確率分布を修正する方法

ということになります。初めに適当に期待値の確率分布を決めて、後から毎回修正を加えていきます。

期待値を修正してくれる尤度関数

ではどのように期待値を修正すれば良いでしょうか？

修正するための材料として持っているのは、アームを引いて得られた結果（報酬）だけです。

報酬を見て、それっぽい期待値に修正する必要があります。実は、この「それっぽい」を表現できる関数があり、尤度関数 と呼ばれています。

尤度関数は、「確率分布の関数形はわかっているが、（期待値、分散などの）母数が不明」というシチュエーションで用いられます。

未知の母数 $\theta$ を持つ密度関数 $ f(x ; \theta)$ *1に対して、観測値 $x_0$ が得られたとき、尤度関数 $L(\theta \mid x_0)$ *2は次のように定義されます。

\begin{aligned} \displaystyle L(\theta \mid x_0) = f(x_0 ; \theta) \end{aligned}

尤度関数は、観測値 $x_0$ が与えられたとき、母数が $\theta$ である尤もらしさ を意味しています。つまり尤度関数の数値が最大となるような $\theta$ は、最も妥当性の高い母数といえます。

トンプソン抽出の文脈に置き換えてみます。

• $\theta_i$ : アーム $i$ の報酬の期待値
• $f(x ; \theta_i)$ : アーム $i$ の報酬が従う確率分布
• $L(\theta_i \mid x_0)$ : アームを引いた結果 $x_0$ に対する現在のアーム $i$ の報酬期待値の尤度関数

とおきます。いま、例えばアーム $i$ の報酬が二項分布に従うことがわかっており、その観測値 $r$ が与えられているとき、尤度関数は、

L(\theta_i \mid r) = {}_n \mathrm{C}_r\theta_{i}^{r} (1 - \theta_i)^{n-r}

となります。広告レコメンドにおけるバンディットアルゴリズムでは、アームの報酬 = 広告のクリック率 として考えることが多いため、今回は二項分布を例に議論していきますが、尤度関数の定義式から明らかなように、尤度関数は報酬の従う確率分布と同形になります。例えば報酬が正規分布に従うならば尤度関数の形も正規分布になりますし、報酬がガンマ分布に従うなら尤度関数もガンマ分布の形になります。

共役事前分布

尤度関数は期待値の妥当性を判定してくれそうなので、期待値の確率分布を修正できそうです。

事前分布を $\pi(\theta)$ 、観測値を $x$、 事後分布を $\pi(\theta \mid x)$ とし、尤度関数を $L(\theta \mid x)$ とすると、ベイズの定理（連続型）から、以下の式が得られます*3:

\begin{aligned} \displaystyle \pi(\theta \mid x) = \frac{\pi(\theta) L(\theta \mid x)}{\int_{\theta}\pi(\theta) L(\theta \mid x)d\theta} \end{aligned}

この辺りはベイズ統計学の基礎的な議論*4なので詳細な議論は割愛しますが、条件付き確率の公式

\begin{aligned} \displaystyle P(A \mid B) = \frac{P(A \mid B) P(B)}{P(A)} \end{aligned}

の派生ぐらいの認識で良いかと思います。

これをトンプソン抽出の文脈に置き換えると、

• $\pi(\theta_i)$ : アームを引く前のアーム $i$ の期待値の分布
• $\pi(\theta_i \mid x_0)$ : アームを引いて $x_0$ を観測した後のアーム $i$ の期待値の分布

となります。先の例で得られた尤度関数 $L(\theta_i \mid r)$ の例だと、アームの期待値 （= クリック率）$\theta_i$ は $0 \leq \theta_i \leq 1$ ですから、アーム $i$ の期待値の事後分布は、

\begin{aligned} \displaystyle \pi(\theta_i \mid r) = \frac{\pi(\theta_i) L(\theta_i \mid r)}{\int_{0}^{1}\pi(\theta_i) L(\theta_i \mid r)d\theta_i} \end{aligned}

となります。

これで期待値の確率分布を修正する手順が揃ったので、あとは「アームの期待値の確率分布」を適当に決めてアルゴリズムをスタートさせたいと思います。が、実際には適当に決めるのではなく、尤度関数に対応する共役事前分布とするのが望ましいとされています。 共役事前分布とは、パラメータを持つ分布で、尤度関数をかけて事後分布を計算すると、（パラメータを除く）関数形が事前分布と一致するような確率分布のことです。

例えば先の例で得られた二項分布の尤度関数

L(\theta_i \mid r) = {}_n \mathrm{C}_r\theta_{i}^{r} (1 - \theta_i)^{n-r}

に対する共役事前分布を考えます。まず、ベイズの定理より、

\begin{aligned} \pi(\theta_i \mid r) &= \frac{\pi(\theta_i) {}_n \mathrm{C}_r\theta_{i}^{r} (1 - \theta_i)^{n-r} }{\int_{0}^{1}\pi(\theta_i) {}_n \mathrm{C}_r\theta_{i}^{r} (1 - \theta_i)^{n-r} d\theta_i} \\ &= \frac{\pi(\theta_i) \theta_{i}^{r} (1 - \theta_i)^{n-r} }{\int_{0}^{1}\pi(\theta_i) \theta_{i}^{r} (1 - \theta_i)^{n-r} d\theta_i} \\ \end{aligned}

ここで、天下り的ですが、変数 $x$ $(0 \leq x \leq 1)$ と2つのパラメータ $a \gt 0, b \gt 0$ を持つような関数

\begin{aligned} f(x; a, b) = \frac{x^{a-1} (1 - x)^{b-1}}{\int_{0}^{1}x^{a-1} (1 - x)^{b-1}dx} \\ \end{aligned}

を考えます。分母を

\begin{aligned} B(a, b) = \int_{0}^{1}x^{a-1} (1 - x)^{b-1}dx \\ \end{aligned}

とおくと、$B(a, b)$ は $x$ に関して定数であるので、$f(x; a, b)$ の全区間 $0 \leq x \leq 1$ での積分値は $1$ になります。また、区間 $0 \leq x \leq 1$ で $f(x; a, b)$ は常に正なので、$0 \leq s \leq t \leq 1$ となる任意の区間 $s \leq x \leq t$ で $f(x; a, b)$ の積分値は正になります。つまり、$f(x; a, b)$ は確率密度関数としての性質を満たします。また、

\begin{aligned} f(x; a+r, b+n-r) &= \frac{x^{a+r-1} (1 - x)^{b+n-r-1}}{\int_{0}^{1}x^{a+r-1} (1 - x)^{b+n-r-1}dx} \\ &= \frac{x^{a-1} (1 - x)^{b-1}x^{r}(1-x)^{n-r}}{\int_{0}^{1}x^{a-1} (1 - x)^{b-1}x^{r}(1-x)^{n-r}dx} \\ &= \frac{\frac{x^{a-1} (1 - x)^{b-1}}{B(a, b)}x^{r}(1-x)^{n-r}}{\int_{0}^{1}\frac{x^{a-1} (1 - x)^{b-1}}{B(a, b)}x^{r}(1-x)^{n-r}dx} \\ &= \frac{f(x; a, b)x^{r}(1-x)^{n-r}}{\int_{0}^{1}f(x; a, b)x^{r}(1-x)^{n-r}dx} \\ \end{aligned}

であるので、先の式に $\pi(\theta_i) = f(\theta_i; a, b)$ を代入すると、$\pi(\theta_i \mid r) = f(\theta_i; a+r, b+n-r)$ が得られます。これより、$f(x; a, b)$ は二項分布に対する共役事前分布であることがわかります。

この $f(x; a, b)$ はベータ分布と呼ばれており、$B(a, b)$ はベータ関数と呼ばれています。

なぜ事前分布は共役でなければならないのか

事前分布と事後分布の関数形が一致すると、アルゴリズムが簡単になるというメリットがあります。

初めに戻ると、トンプソン抽出のアルゴリズムは、

1. それぞれのアームの期待値の確率分布に従って乱数を生成する
2. 生成された乱数が最大のアームを1回引いて期待値の確率分布を修正する

という手順を繰り返すものでした。いま、アーム $i$ $(i = 1, \dots, K)$ の期待値の確率分布に対応する共役事前分布を $f(\theta_i)$ とおきます。ここまでの議論を踏まえると、これらは具体的に次のような手順に落とすことができます:

1. 各$i$ $(i = 1, \dots, K)$ について、$f(\theta_i)$ に従って乱数を生成する
2. 生成された乱数が最大のアーム$i^{*}$を1回引いて事後分布を求め、それを新たに $f(\theta_ {i^{*}})$ とする

アルゴリズム中、前の試行で得られた事後分布は次の試行での事前分布として利用されますが、これらの関数形が異なる場合、アルゴリズムのループの中で毎回関数形を選択し直さなければなりません。共役事前分布の場合、尤度関数をかけても関数形は変わらないので、常に同じ関数を（パラメータだけ変えて）使い回すことができます。そのため、事前分布として初期値を設定する際には、尤度関数に対応する共役事前分布を設定するのが望ましいとされています。

今回の例では、尤度関数が二項分布の形であったため、対応する共役事前分布はベータ分布でしたが、二項分布以外の尤度関数の形にも対応する分布がいくつか知られています。有名なものを以下に挙げておきます。

尤度関数の形	共役事前分布
ベルヌーイ分布	ベータ分布
二項分布	ベータ分布
正規分布(分散が既知)	正規分布
正規分布(分散が未知)	Scaled inverse chi-squared分布
ポアソン分布	ガンマ分布

まとめ

トンプソン抽出の実体は、ベイズの定理による期待値の推定でした。また、共役事前分布を初期値として選ぶことで、アルゴリズムがシンプルな形で実装できるということも確認できました。

トンプソン抽出についての記事や書籍の解説を読んでいて、省略されがちな議論を自分なりに補完したいと思い、今回の記事作成に至りました。最後に練習問題も載せておきます。本記事と併せて、理解の一助となれば幸いです🐟 *5

練習問題

Scaled inverse chi-squared分布とは次のような分布である:

f(x; \nu, \tau^2) = \left(\frac{\nu\tau^2}{2}\right)^{\frac{\nu}{2}}\frac{e^{-\frac{\nu\tau^2}{2x}}}{\Gamma(\frac{\nu}{2})x^{1+\frac{\nu}{2}}} \hspace{10pt}(\nu > 0, \tau > 0)

ただし、$\Gamma(z)$ はガンマ関数であり、$z > 0$ に対して次のように定義されている。

\Gamma(z) = \int_{0}^{\infty} t^{z-1} e^{-t}dt

これについて、次の問いに答えよ。

(1) 期待値 $\mu$ が未知であり、分散 $\sigma^{2}$ が既知の正規分布 $N(x; \mu, \sigma^{2})$ に従うアームを1回ひき、報酬 $r$ を観測した。このとき、 $\mu$ に関する尤度関数 $L(\mu \mid r)$ を求め、正規分布 $N(\mu; \eta, \tau^{2})$ が共役事前分布となることを示せ。

(2) 期待値 $\mu$ が既知であり、分散 $\sigma^{2}$ が未知の正規分布$N(x; \mu, \sigma^{2})$に従うアームを1回ひき、報酬 $r$ を観測した。このとき、$\sigma^{2}$ に関する尤度関数 $L(\sigma^{2} \mid r)$ を求め、Scaled inverse chi-squared分布 $f(\sigma^{2}; \nu, \tau^{2})$ が共役事前分布となることを示せ。

(3) n個の観測値 $x_1, \dots, x_n$ に対して、尤度関数は一般に次のように定義される:

\begin{aligned} \displaystyle L(\theta \mid x_1, \dots, x_n) = \prod_{i = 1}^{n}f(x_i ; \theta) \end{aligned}

　いま、(2)の正規分布に従うアームを $n$ 回ひき、報酬 $r_1, \dots, r_n$ を観測した。このとき、 $\sigma^{2}$ に関する尤度関数 $L(\sigma^{2} \mid r_1, \dots, r_n)$ を求め、Scaled inverse chi-squared分布 $f(\sigma^{2}; \nu, \tau^{2})$ が共役事前分布となることを示せ。

最後に

Speeeでは一緒にサービス開発を推進してくれる仲間を大募集しています！
こちらのFormよりカジュアル面談も気軽にお申し込みいただけます！

Speeeでは様々なポジションで募集中なので「どんなポジションがあるの？」と気になってくれてた方は、こちらチェックしてみてください！もちろんオープンポジション的に上記に限らず積極採用中です！！！

*1:$f(x ; \theta)$は「変数を $x$ とし、定数のパラメータ$\theta$を持つ関数 $f$」という意味で使っています。

*2:$L(\theta \mid x_0)$ は、条件付き確率の記法にならい、「観測値 $x_0$が得られたときの、母数 $\theta$ の尤度関数 $L$ 」という意味で使っています。尤度関数は確率分布ではないので、この記法は記号の濫用ですが、直感的であり広く採用されている記法なので、本記事でもこの記法に従います。

*3:複数の観測値を用いてベイズの定理から事後分布を計算するには、それぞれの観測値が i.i.d. （独立かつ同一分布）であるという条件が必要です。このようなケースについては、最後の練習問題 (3) で触れています。

*4:余談ですが、一般的な数理統計学の本には、ベイズの定理、尤度関数、最尤推定法までは載っていても、ベイズの定理の連続型の公式は載っていないことが多く、私は「ベイズ」統計学の本を読んで初めて出会いました。ベイズの定理の連続型の公式は尤度関数の強力な応用ですし、議論としても地続きに学習できる内容だと思います。それなのに一般的な数理統計学の本で扱われていないのはなぜなのでしょうか…

*5:本文の修正および練習問題解答などのリクエストは、Twitter @mov_vc で連絡いただけると助かります。