備忘録

主に機械学習、統計、競技プログラミングに関する備忘録

Robust Bayesian inference via coarsening

Robust Bayesian inference via coarsening についてまとめる。

背景

通常のベイズ推定ではあるデータの生成過程を仮定することで推論を行うが、これらの仮定が実際の生成過程と差異がある場合、僅かな差異であっても推定結果に大きな違いをもたらすことが知られている。
この論文ではこういった差異に対してロバストな結果を与えることが出来る coarsened posterior (または c-posterior)というものを提案している。

Toy Example

$X_1, \cdots, X_n$ がそれぞれ $\mathrm{Bernoulli}(\theta)$ からの i.i.d. サンプルであるとする。次に $H_{0} : \theta=1/2$ versus $H_{1} : \theta \ne 1/2$ の検定を考える。
各仮設に対する事前分布は一様に $\Pi(H_0) = \Pi(H_1) = 1/2$ であり、$H_1$ の下での $\theta$ の事前分布は $\theta \mid H_1 \sim \mathrm{Uniform}(0, 1)$ に従うものとする。
このとき、観測データ $x_{1:n}$ が実際に $\mathrm{Bernoulli}(\theta)$ からのサンプルであれば、事後分布は正しい "答え" を導くことが保証されてる。すなわち、

$$ \Pi(H_0 \mid x_{1:n}) \overset{\mathrm{a.s.}}{\rightarrow} \mathbb{1}(\theta = 1/2) $$

が成り立つ。ここで$\mathbb{1}(\cdot)$ は $E$ が真であれば $\mathbb{1}(E) = 1$、そうでないとき $\mathbb{1}(E) = 0$ となるような指示関数である。
しかし、現実的な問題として、仮定した生成過程が厳密に正しいケースは稀である。例えば $\theta = 0.51, 0.56$ と設定して事後分布を見てみると、従来の事後分布では $n$ が大きくなるにつれ所望の値から大きく乖離することが確認できる。

f:id:invariant:20200815211407p:plain

上の例の場合 $\theta = 0.5$ と $\theta=0.51$ の差も検出したいのであれば問題はないのだが、あるモデルを仮定したもとで予測を行いたい場合、真のモデルと僅かに乖離しただけで予測精度が大幅に悪化するのは避けたいため、このように使用するモデルに対してロバストな事後分布を構成することには大きなご利益がある。

coarsened posterior

本論文の筆者はモデルの設計者が決めたパラメトリックモデルにおける確率分布 $P_{\theta_I}$ と実際に観測されるデータの経験分布 $\hat{P}_{x_{1:n}}$ の間の距離(ダイバージェンス)に関する条件付け行う。具体的な定義を行うため、記号の準備を行う。
理想的なモデルを生成する分布(idealized distribution)を含むパラメトリックモデルを $\left\{P_\theta : \theta \in \Theta\right\}$ とする。特に、idealized distributionを実現するパラメータを $\theta_I$ とおく。$P_{\theta_I}$ からの現実には観測されないサンプルを $X_1, \cdots, X_n \in \mathcal{X}$ とする。次に、実際に観測されるサンプルを $x_1, \cdots, x_n \in \mathcal{X}$、経験分布を $\hat{P}_{x_{1:n}} = \frac{1}{n} \sum_{i=1}^{n} \delta_{x_{i}}$ とおく。最後に、2つの確率分布の距離を測る関数として $d(\cdot, \cdot)$ を用意する。通常のベイズ推定では $X_{1:n} = x_{1:n}$ であると仮定することを確認しておく。
上でも述べたように、一般に観測値が設計者の設定した分布に正確に従うことはありえないので、$X_{1:n}$ の仮想的な経験分布 $\hat{P}_{X_{1:n}}$ と実際の観測値の経験分布 $\hat{P}_{x_{1:n}}$ 間の距離に一定の"遊び"が入るのが自然と考えられる。
例えば両確率分布間に距離 $r$ までを許容すると決めたならば、条件 $d_n(X_{1:n}, x_{1:n}) := d(\hat{P}_{X_{1:n}}, \hat{P}_{x_{1:n}}) < r$ の下での $\theta$ の事後分布を構成するのが良さそうである。また、$r$の具体的な値を決めることは困難なので、確率変数化し、 $R \sim H$ に従うものとする。
以上をまとめるとc-posteriorは以下のように定義される。

Definition( c-posterior ) $$ \begin{align} \Pi(d\theta \mid d_n(X_{1:n}, x_{1:n})) &\propto \Pi(d\theta) \mathbb{P}(d_n(X_{1:n}, x_{1:n}) < R \mid \theta) \\ &=: \Pi(d\theta) \int_{\mathcal{X}^n} G(d_{n}(x'_{1:n}, x_{1:n})) P_\theta^{n}(dx'_{1:n}) \end{align} $$

ここで、$G(\cdot)$ は $G(r) = \mathbb{P}(R > r)$ を満たす関数である。 この積分はまず $x'_{1:n}$ が $P^{n}_\theta$ からサンプルされ、次に確率変数 $R$ が $d_n$ を上回る確率を畳み込んだものと読むとよい。これはある $\theta$ の下で観測値 $x_{1:n}$ とそれなりに近いサンプルがどれだけ生成され得るかを表すような値と解釈出来る。
このままでは具体性にかけるので、 $d, H$ を具体的に定めて計算を進める。$d$ をKLダイバージェンス、$R \sim \mathrm{Exp}(\alpha) \propto \exp(-\alpha R)$ とすると、c-posteriorは $G$ に指数分布の密度関数を代入することで、

$$ \begin{align} \Pi(d\theta \mid d_{n}(X_{1:n}, x_{1:n}) < R) &\propto \mathbb{E}\left[ \exp(-\alpha D(\hat{p}_{x_{1:n}} \mid\mid \hat{p}_{X_{1:n}} ) ) \mid \theta\right] \Pi(d\theta) \\ &\propto \exp(\zeta_n \sum_{i=1}^{n} \log p_{\theta} (x_i) ) \Pi(d\theta) \\ &= \Pi(d\theta) \prod_{i=1}^{n} p_{\theta}(x_i)^{\zeta_n}, \\ \zeta_n &= \frac{1}{1 + n / \alpha} \end{align} $$

が得られる。ただし、第1行から第2行へはそれなりの飛躍あるため、詳細はAppendixを参照。ここまで来れば通常の尤度に逆温度 $1/\xi_i$ を導入した事後分布であると解釈でき、MCMCなどの各種計算も可能であるが分かる。既存の数値計算アルゴリズムをそのまま流用してロバストな事後分布が構成出来るというのがこの論文の貢献である。

理論的結果

事後分布が真の生成分布 $P_{o}$ の変化に対してロバストであることは、事後分布が(何らかの位相において)$P_o$ に対して連続的であると言える。これはラフに言えば $P_o$ の微小な変化に対して事後分布も微小な変化をすることを要請する。
論文中では以下の定理で c-posterior が $P_o$ に関して連続であることを示している。

回帰への応用

論文中では時系列モデルと回帰モデルへの応用を紹介している。ここでは後者の結果について紹介しよう。
モデルとしてはSpike-slab modelによるスパースな回帰モデルを考える。

$$ \begin{align} W &\sim \mathrm{Beta}(r, s) \\ \beta_{j} &\sim \mathcal{N}(0, 1/L_0) \mathrm{with probability} W, \mathrm{otherwise} \beta_j = 0, j=1,\cdots,p \\ \lambda &\sim \mathrm{Gamma}(a, b)\\ Y_i \mid \beta, \lambda &\sim \mathcal{N}(\beta^\top x, 1/\lambda), i=0,\cdots,n \end{align} $$

これに対して、以下の分布に従うデータ $(x_i, y_i)$ を生成させる。

$$ y_i = \beta_{01} + \beta_{02}(x_{i2} + \frac{1}{16} x_{i2}^2) + \epsilon_{i} $$

ただし、 $\beta_{01} = -1, \beta_{02} = 4, \epsilon \sim \mathcal{N}(0, 1)$ とする。

f:id:invariant:20200816190015p:plain

このデータにして事後分布 $\pi(\beta, \lambda \mid y) \propto \pi(\beta, \lambda) p(y \mid \beta, \lambda)^{\zeta_{n}}$ を計算し、非ゼロ成分の個数と $\beta$ の各成分の累積分布関数を求めると以下のようになる。

f:id:invariant:20200816190500p:plain

まず非ゼロ成分の個数について、通常のベイズ事後分布ではサンプルサイズが大きくなるにつれ3~5と予測するようになるが、c-posteriorではサンプルサイズが増えても正しい非ゼロ成分の個数を当てていることが分かる。
さらに、事後分布の最頻値を見ると、$\beta_{1}$ が$-0.7$ 付近、$\beta_2$ が $3.5 \sim 4.0$ の間と、実際の値に近い値を取ってることが分かる。

まとめと感想

本論文で提案された事後分布は非常にシンプルな形で、既存のプログラムにも組み込みやすそうである。 数値実験の結果についても、真の分布がわからない場合により妥当な推論が可能になることが確認出来た。時間があれば実際に手元でデータを動かしてみて、その有用性を確認してみたい。
一点気になることとしては、回帰などの予測が伴うタスクの際に、c-posteriorを使うことでどの程度予測精度が改善されるかが分からないので、調べてみたい。