Randomized Response のベイズ推論(3): 変分ベイズ

Randomized Response はアンケートの回答をランダム化することで、個人の回答は伏せつつ平均などの統計量を得る手法の1つ。
前回までの記事で、ランダム化された回答から真の割合を最尤推定ベイズ推定(ギブスサンプリング)で推定する方法とそれらの実験を紹介した。今回は同じベイズ推定でも変分ベイズによる推論を行ってみる。
Randomized Response のギブスサンプリングによる推論までは普通に論文が見つかるが、変分ベイズ推論は探した範囲では見つけられなかったので、頑張って導出してみる。

Collapsed Variational Bayes 推論

記法は前回と同様。
D 通りの選択肢を持つ質問 Q に対する回答を確率変数 X\in\{1,\cdots,D\} で表す。
回答者は真の回答 X を知られたくないので、 X をランダム化した Y\in\{1,\cdots,D\} を集計者に返す。
このとき以下のモデルを考える。

 \begin{cases}
P(Y=j|X=i)=p_{ij} & i,j=1,\cdots,D, \; \sum_j p_{ij}=1 \\
P(X=i)=\pi_i & i=1,\cdots,D \\
P(\boldsymbol{\pi}) = \rm{Dir}( \boldsymbol{\pi}| \boldsymbol{\alpha} )
\end{cases}

ランダム化の確率 P(Y=j|X=i)=p_{ij} からなる計画行列を \boldsymbol{P}=\left(p_{ij}\right) と記す。
また \boldsymbol{\pi}=(\pi_i) にはパラメータ \boldsymbol{\alpha}=\left( \alpha_i \right) を持つディリクレ共役事前分布を入れる。

ユーザ n の真の回答を X_n、それをランダム化したものを Y_n で表す。
集計者は N 人の回答者のランダム化された回答の観測値 Y_1,\cdots,Y_N から \mathbb{E}[X] などを推定する。

f:id:n_shuyo:20210119180208p:plain:w320

変分ベイズ推論は事後分布に独立性の仮定(変分近似)を入れて、最終的に \log q(X_n)\log q(\boldsymbol{\pi}) などを計算して更新式を導出する

 q(X_1,\cdots,X_N,\boldsymbol{\pi}) \approx q(\boldsymbol{\pi})\prod_{n=1}^Nq(X_n)

しかしせっかくなので、より高い性能を期待して LDA(Latent Dirichlet Allocation) でおなじみの Collapsed Variational Bayes を使ってみよう(Teh+2007, Asuncion+2009)。
Collapsed Variational Bayes ではまず変分近似を以下のように弱める。

  q(\boldsymbol{X},\boldsymbol{\pi}) 
= q(\boldsymbol{\pi}|\boldsymbol{X})q(\boldsymbol{X})
\approx q(\boldsymbol{\pi}|\boldsymbol{X})\prod_{n=1}^Nq(X_n)

ただしこれ以降 X_1,\cdots,X_N\boldsymbol{X}Y_1,\cdots,Y_N\boldsymbol{Y}X_n 以外の X_* 全体を \boldsymbol{X}^{-n}=\{X_1,\cdots,X_{n-1},X_{n+1},\cdots,X_N\} で表す。
この仮定に基づき変分自由エネルギー \mathcal{F}\left(q(\boldsymbol{X},\boldsymbol{\pi})\right) を以下のように展開する。

 \begin{eqnarray}
\mathcal{F}\left(q(\boldsymbol{X},\boldsymbol{\pi})\right)
&=& \mathbb{E}_{q(\boldsymbol{X},\boldsymbol{\pi})} [ -\log P(\boldsymbol{Y},\boldsymbol{X},\boldsymbol{\pi}) ] - \mathcal{H}(q(\boldsymbol{X},\boldsymbol{\pi})) \\
&=& \mathbb{E}_{q(\boldsymbol{X})} [ \mathbb{E}_{q(\boldsymbol{\pi}|\boldsymbol{X})} [ -\log P(\boldsymbol{Y},\boldsymbol{X},\boldsymbol{\pi}) ] - \mathcal{H}(q(\boldsymbol{\pi}|\boldsymbol{X})) ]
 - \mathcal{H}(q(\boldsymbol{X}))
\end{eqnarray}

ここで自由エネルギー \mathcal{F}\left(q(\boldsymbol{X},\boldsymbol{\pi})\right)q(\boldsymbol{X},\boldsymbol{\pi}) で最小化する代わりに、まず q(\boldsymbol{\pi}|\boldsymbol{X}) で最小化して、次に q(\boldsymbol{X}) で最小化する。q(\boldsymbol{\pi}|\boldsymbol{X}) で最小化したものを \mathcal{F}\left(q(\boldsymbol{X})\right) と書く。

 \begin{eqnarray}
\mathcal{F}\left(q(\boldsymbol{X})\right) 
&:=& \min_{q(\boldsymbol{\pi}|\boldsymbol{X})} \mathcal{F}\left(q(\boldsymbol{X},\boldsymbol{\pi})\right) \\
&=& \mathbb{E}_{q(\boldsymbol{X})} [  -\log P(\boldsymbol{Y},\boldsymbol{X}) ]
 - \mathcal{H}(q(\boldsymbol{X}))
\end{eqnarray}

この \mathcal{F}\left(q(\boldsymbol{X})\right) を、q(\boldsymbol{X})=\prod_{n=1}^Nq(X_n) に基づいて、各 q(X_n) で反復的に最小化するのが Collapsed Variational Bayes となる。
多項分布となる q(X_n) のパラメータを q(X_n=i)=\gamma_{ni} とし、c_k^{-n}=\sum_{m\neq n}I(X_m=k), \;\gamma_k^{-n}=\sum_{m\neq n} \gamma_{mk} とおくと(ただし I(\cdot) は指示関数)、その更新式は以下のように求められる。

 \begin{eqnarray}
\gamma_{ni} &\propto& \exp\left( \mathbb{E}_{q(\boldsymbol{X}^{-n})} [ \log P(\boldsymbol{Y},\boldsymbol{X}^{-n},X_n=i) ]\right)  \\
&\propto& \exp\left(  \log P(Y_n|X_n=i) + \mathbb{E}_{q(\boldsymbol{X}^{-n})} \left[ \log \int \pi_i \prod_k \pi_k^{\alpha_k+c_k^{-n}-1}d\boldsymbol{\pi}  \right]\right)  \\
&\propto& \exp\left(  \log p_{iY_n} + \mathbb{E}_{q(\boldsymbol{X}^{-n})} \left[ \log (\alpha_i+c_i^{-n})  \right]\right)  \\
&\approx& p_{iY_n}(\alpha_i+\gamma_i^{-n})
\end{eqnarray}

式の最後の近似は、LDA の CVB0 と同じくテイラー展開の 0 次の項のみに近似している。

実験

実装はこちら。

github.com

同じ反復法だが、データ件数 N でループする必要があるため、Gibbs Sampling より処理が10倍以上重い。そのため、推論の試行とサマリーを分離した上、試行の並列化までしている(苦笑)。
ここまでやっても、(N,\alpha) の組み合わせ 9通り×10000試行の推論には手元の環境で10時間かかった。

前回と同様に Randomized Response の計画行列 P(Y|X) は1個のパラメータ p で表現し、実験では p=1/5 を用いる。
 \displaystyle P(Y=j|X=i)=p_{ij}=\begin{cases}p + \frac1D(1-p)& (i=j) \\ \frac1D(1-p) & (i\neq j) \end{cases}

D=4 個の選択肢を持つ質問に対し、N=100,1000,10000 人のユーザが 1:2:3:4 の割合で回答(\boldsymbol{\pi}=(0.1\;0.2\;0.3\;0.4)^\top)p=1/5 の Randomized Response でランダム化された回答から推定した \hat\pi_1,\cdots,\hat\pi_4 の分布を見る。
事前分布のパラメータによる挙動の違いも見たいので、\alpha=1.0,\;0.1,\;0.01 を対称なディリクレ事前分布のパラメータとして用いる。

まずは手法ごとの違いを見るため、最尤推定、ギブスサンプリング、そして変分ベイズ(Collapsed Variational Bayes)に対して N=100,1000,10000 ごとの推定値の分布を見比べてみる。事前分布のパラメータは \alpha=1 のものを選んでいる。
今回は密度推定した線も書き加えてみたので、分布の形がよりわかりやすい。
 

N=100 N=1000 N=10000
MLE f:id:n_shuyo:20210204133030p:plain:w200 f:id:n_shuyo:20210204133034p:plain:w200 f:id:n_shuyo:20210204133037p:plain:w200
Gibbs f:id:n_shuyo:20210204130853p:plain:w200 f:id:n_shuyo:20210204130856p:plain:w200 f:id:n_shuyo:20210204130900p:plain:w200
VB f:id:n_shuyo:20210204130903p:plain:w200 f:id:n_shuyo:20210204130907p:plain:w200 f:id:n_shuyo:20210204130911p:plain:w200

最尤推定解は \pi_i によらず均等な散らばり具合で [0,1] からはみ出すが、ベイズ化によって \pi_i を確率値としてモデリングすると解消できることは前回までの記事で見たとおり。変分ベイズももちろんベイズモデルなので、同様に [0,1] に収まるように推定される。
上の図ではギブスサンプリングと変分ベイズの結果はそっくり同じに見えるが、実は縦軸・横軸の数値をよく見ると分かるように、変分ベイズの方が分散が小さい(山も高い)。
それがもっと明確に読み取れるように、\hat\pi_1\hat\pi_4 について箱ひげ図を書いて、特にそれらの散らばり具合の違いを見比べてみる。

N=100 N=1000 N=10000
f:id:n_shuyo:20210204152557p:plain:w200 f:id:n_shuyo:20210204152600p:plain:w200 f:id:n_shuyo:20210204152603p:plain:w200

すべての場合で 変分ベイズ(左) のほうがギブスサンプリング(右) より分散が小さい。横軸の下に書かれた推定値の標準偏差からも、グラフの見た目だけではなく実際に分散が小さくなっていることがわかる。これは推定値の精度が良くなっていることを意味している。*1

次にディリクレ事前分布のハイパーパラメータによる挙動の違いを見る。
N=10000 ではハイパーパラメータによる差異がほとんどない(データが多いと事前分布の影響が小さくなる)のは前回ギブスサンプリングでも見たとおりなので、N=100,1000 の図を掲載する。

N=100 α=0.01 α=0.1 α=1.0
Gibbs f:id:n_shuyo:20210204161937p:plain:w200 f:id:n_shuyo:20210204161941p:plain:w200 f:id:n_shuyo:20210204161944p:plain:w200
VB f:id:n_shuyo:20210204161959p:plain:w200 f:id:n_shuyo:20210204162002p:plain:w200 f:id:n_shuyo:20210204162005p:plain:w200
N=1000 α=0.01 α=0.1 α=1.0
Gibbs f:id:n_shuyo:20210204162023p:plain:w200 f:id:n_shuyo:20210204162026p:plain:w200 f:id:n_shuyo:20210204162030p:plain:w200
VB f:id:n_shuyo:20210204162047p:plain:w200 f:id:n_shuyo:20210204162051p:plain:w200 f:id:n_shuyo:20210204162054p:plain:w200

ギブスサンプリングでは \alpha を 1 より小さくすると事前分布の影響が強く出た推定になっていた。特に N=100,\alpha=0.01 の両端が上がった分布(つまり \hat\pi_i のどれか1つが 1 で残りが 0)が顕著だろう。*2
それに対し、変分ベイズではそうした極端さはほぼなく、直感に即した真値の周りの山形の分布となっている。

\alpha を小さくすることによる平滑化の低減(トレードオフとして分散の増加)を箱ひげ図で確認しよう。平均が理屈通りに振る舞うのは当然なので、今回は中央値に着目してみる。

α=0.01 α=0.1 α=1.0
N=100 f:id:n_shuyo:20210204165425p:plain:w200 f:id:n_shuyo:20210204165428p:plain:w200 f:id:n_shuyo:20210204152557p:plain:w200
N=1000 f:id:n_shuyo:20210204165531p:plain:w200 f:id:n_shuyo:20210204165534p:plain:w200 f:id:n_shuyo:20210204152600p:plain:w200

まず前回の記事で指摘しそびれていたこととして、ギブスサンプリングの \alpha<1 での極端な分布では中央値も真値から大きく離れていることが箱ひげ図から読み取れる。推定値の分散も非常に大きい。N=100,\alpha=0.01\hat\pi_4 では台が [0,1] しかないのに標準偏差が 0.430 と、一様分布(標準偏差 0.289)のほうがマシなレベル。
ギブスサンプリングで \alpha<1.0 のハイパーパラメータを使うのは難しそうだ。

変分ベイズでは、真値の周りの山型分布から想定できるとおり、中央値がちゃんと真値に近い値として推定できている。そして \alpha を小さくすると中央値はより真値に近づく。

このように変分ベイズ(Collapsed Variational Bayes)による推論は、最尤推定やギブスサンプリングより良い推定値が得られる。\alpha<1 でも直感的に良い推定値が得られる可能性が十分高く、欲しい精度に合わせてハイパーパラメータを厳選する甲斐がありそう。
実際の観測値に対する推定を行う場合、平均より中央値が外れてないことのほうが嬉しいので(確率 1/2 で中央値以上/以下になる)、その点でも変分ベイズによる推論は筋が良いとわかる。

同じ(真の事後)分布の推定を目的としているギブスサンプリングとの意外なほどの大きな差異が何によって生じているのかは興味深いところ。思い当たる要因は変分近似と CVB0(テーラー展開の0次近似)しかなく、どちらも原理的には精度を下げる効果しかないはずなのに、変分ベイズのほうが自明に良い結果を叩き出すのは、LDA-CVB0 の不思議な性能の良さと同じ根っこにつながっていそう(そう言えば CVB0 の理論的解釈とか誰かやってたような……)。

References

  • Teh, Yee W., David Newman, and Max Welling. A collapsed variational Bayesian inference algorithm for latent Dirichlet allocation. In Advances in Neural Information Processing Systems 19, 2007.
  • Asuncion, A., Welling, M., Smyth, P., and Teh, Y. W. On smoothing and inference for topic models. In Proceedings of the International Conference on Uncertainty in Artificial Intelligence, 2009.

*1:ギブスサンプリングはサンプルの平均を推定値とするが、その長さが小さいために散らばりが大きくなっているという可能性を確認するため、元のソースでは 200サンプルを burn-out して、続く 200サンプルの平均をとっていたのを 500サンプル burn-out & 500サンプル平均 にしてみたが、傾向は全く変わらなかった。

*2:それでいて \mathbb{E}[\hat\pi_i] はちゃんと真値に近くなるのがよくできているというかなんというか