Randomized Response のベイズ推論(2): ギブスサンプリング

Randomized Response はアンケートの回答をランダム化することで、個人の回答は伏せつつ平均などの統計量を得る手法の1つ。
前回記事では、回答の割合の推定量最尤推定で得る手順を紹介したが、割合の推定値が負になる可能性があることを示した。
shuyo.hatenablog.com
その問題を解消する方法の1つに、統計モデルをベイズ化するアプローチがある。今回はその中でも一番シンプルなギブスサンプリングによる推論を紹介する。

2値アンケートだとベイズ化の効果が見えにくいので、ベイズ化の前にまずは多値化する。
なお、i 番目の要素が a_i であるベクトルを \left(a_i\right)(i,j)要素が a_{ij} である行列を \left(a_{ij}\right) のように簡易的に表記している。

質問 Q に対する回答者の答えを確率変数 X で表す。
回答者は真の回答 X を知られたくないので、 X をランダム化した Y を集計者に返す。
集計者は N 人の回答者のランダム化された回答 Y_1,\cdots,Y_N から \mathbb{E}[X] などを推定する。

QD 通りの回答を持つとし、X, Y ともに 1,\cdots,D の値を取るとする。
X をランダム化した Y の確率を P(Y=j|X=i)=p_{ij} とし\left(\sum_j p_{ij}=1\right)p_{ij} を要素とする計画行列 \boldsymbol{P}=\left(p_{ij}\right) を考える。
X のモデルを P(X=i)=\pi_i としたとき、\boldsymbol{\pi}=\left(\pi_i\right) の推定値が欲しい回答の割合となる。

Y の観測値 \mathcal{D}=\left\{y_1,\cdots,y_N\right\} に対し、過程は省略するが、最尤推定\hat{\boldsymbol{\pi}} は以下のように求められる。ただし c_j=\left|\left\{ n | y_n=j \right\}\right|, \boldsymbol{c}=\left(c_j\right) とする。

\displaystyle \hat{\boldsymbol{\pi}}=\boldsymbol{P}^{-1}\boldsymbol{c}

計画行列 \boldsymbol{P}=\left(p_{ij}\right)D(D-1) 個のパラメータを持つが、自由度が多すぎてそのままでは使いにくい。
そこで真値が X=i であるとき、確率 pY=i、確率 1-p\{1,\cdots,D\} から一様乱数で選んだ値を Y とするランダム化メカニズムがよく採用される。
この場合、1つのパラメータ p に対し、計画行列は以下のように定められる。*1

\displaystyle p_{ij}=\begin{cases}p + \frac1D(1-p)& (i=j) \\ \frac1D(1-p) & (i\neq j) \end{cases} (式1)

最尤推定\hat{\boldsymbol{\pi}} は 0 から 1 の範囲に収まらないことがある。前回と同様にシミュレーションで確認しておこう。
実験に使ったコードはこちら:
github.com

4個の選択肢を持つ質問に対し、N(=100,1000,10000) 人のユーザが 1:2:3:4 の割合で回答するとき、ランダム化された回答から \hat{\boldsymbol{\pi}} を推定する。
p=1/5 の Randomized Response でランダム化された回答から推定した \hat\pi_1,\cdots,\hat\pi_4 をプロットした。

f:id:n_shuyo:20210119133717p:plain:w320f:id:n_shuyo:20210119133720p:plain:w320
f:id:n_shuyo:20210119133723p:plain:w320

N=10000 ではいい感じに推定できているが、N=100 では混然一体となってしまっており、\pi_i 間の大小さえうまく推定できるかどうかわからないことが読み取れる。
Randomized Response のパラメータ p を大きくすれば精度は改善できるので、集めたい N の大きさや最低限望む安全性などと相談して p を決める必要がある。

N=100 における各 \pi_i の平均・標準偏差、および 95%信頼区間も見ておこう。

\pi_i 0.1 0.2 0.3 0.4
\mathbb{E}[\hat\pi_i] 0.1 0.2 0.3 0.4
\rm{std}[\hat\pi_i] 0.21 0.21 0.22 0.22
95%区間 -0.3~0.5 -0.2~0.6 -0.1~0.75 0.0~0.85

これほど重なっていても、\mathbb{E}[\hat\pi_i]=\pi_i はちゃんと成立している。
分散(標準偏差)は i によらず一定だ。ここはベイズ版との比較ポイントになるので覚えておこう。
95%信頼区間は思いっきり負の領域にはみだしている。ベイズ化によってこれの改善を図る。

ベイズ版 Randomized Response

\boldsymbol{\pi} に共役事前分布としてディリクレ分布を入れて、その事後分布を推定する。

\displaystyle P(\boldsymbol{\pi}) = \rm{Dir}( \boldsymbol{\pi}| \boldsymbol{\alpha} ) = \frac{\Gamma\left(\sum_{i=1}^D\alpha_i\right)}{\prod_{i=1}^D\Gamma(\alpha_i)}\prod_{i=1}^D\pi_i^{\alpha_i-1}

\boldsymbol{\alpha}=\left( \alpha_i \right) はディリクレ事前分布のパラメータである。

事後分布の推定には方法はいくつかあるが、ここではまずギブスサンプリングを使った推論を導出しよう。
ユーザ n の真の回答を X_n、それをランダム化したものを Y_n という確率変数で表す。
また X_n 以外の X_* 全体を \boldsymbol{X}_{-n}=\{X_1,\cdots,X_{n-1},X_{n+1},\cdots,X_N\} で表す。

f:id:n_shuyo:20210119180208p:plain:w490
Randomized Response のグラフィカルモデル

このとき観測値 \mathcal{D}=\left\{Y_1=y_1,\cdots,Y_N=y_N\right\} に対する事後分布 P(\boldsymbol{\pi}|\mathcal{D}) をギブスサンプリングで求めるのに必要な、潜在変数 X_1, \cdots, X_N および \boldsymbol{\pi} の全条件付き事後分布を計算する。

\displaystyle \begin{eqnarray}
&& P(X_n=i|\boldsymbol{X}_{-n},\boldsymbol{\pi},\mathcal{D}) \\
&\propto& P(X_n=i,Y_n=y_n|\boldsymbol{\pi}) \\
&=& P(X_n=i|\boldsymbol{\pi}) P(Y_n=y_n|X_n=i) = \pi_i p_{iy_n}
\end{eqnarray} (式2)

\displaystyle \begin{eqnarray}
&& P(\boldsymbol{\pi}|X_1,\cdots,X_N,\mathcal{D}) \\
&\propto& P(\boldsymbol{\pi},X_1,\cdots,X_N,\mathcal{D}) \\
&=& P(\boldsymbol{\pi})\prod_{n=1}^N \left\{ P(X_n|\boldsymbol{\pi})P(Y_n|X_n) \right\}
\propto \prod_{i=1}^D \pi_i^{\alpha_i-1+d_i}
\end{eqnarray} (式3)

ただし d_iX_n の値が i である個数 d_i = \left|\{n|X_n=i\}\right| とする。
P(\boldsymbol{\pi}|X_1,\cdots,X_N,\mathcal{D}) はその式の形から  \rm{Dir}\left( \boldsymbol{\pi}| (\alpha_i+d_i) \right) とわかる。

X_1, \cdots, X_N\boldsymbol{\pi} の初期値を適当に決めて、(式2) と (式3) からのサンプリングを繰り返せばギブスサンプリングとなる。
(式2) は、X_n のサンプリングが他の \boldsymbol{X}_{-n} に依存しないことを表している。実際グラフィカルモデルを見ると、\boldsymbol{\pi} を縛れば X_1,\cdots,X_N は条件付き独立になる。
つまり X_1,\cdots,X_N 全体を一度にサンプリングするようなコードを書くことができる*2

実験に使ったコードはこちら:
github.com

実装ではサンプリングを400回繰り返し、最初の200回を捨てて、残る 200回のサンプル \boldsymbol{\pi} の系列の平均をとって割合の推定値としている(回数は適当w)。
事前分布のパラメータは {}^\forall \alpha_i=\alpha と対称とし、まずは \alpha=1 に対してデータの件数 N=100,1000,10000 に対して、Randomized Response の試行とギブスサンプリングによる推論を繰り返して得られた分布が以下のようになった。

f:id:n_shuyo:20210120155525p:plain:w240f:id:n_shuyo:20210120155534p:plain:w240f:id:n_shuyo:20210120155546p:plain:w240

まず N=10000 ではきれいな結果になっている。これは問題ないだろう。

N が小さくなると分散が大きくなり、最尤推定解が負になることがある問題を解決するためにベイズ化したわけだが、N=100 や 1000 のヒストグラムの横軸を見れば分かる通り、\hat\pi_1 などの分布もすべて 0 以上の範囲に収まっている。
わかりやすい N=1000 で見ると、\hat\pi_2,\hat\pi_3,\hat\pi_4 は左右対称な釣鐘型であるのに対し、\hat\pi_1 の分布は左の裾が 0 のところにストンと落ちる非対称な分布となっている。
N=100 になると分布の重なりは著しくなってしまうが、それでも分布の台が 0 から 1 の範囲に収まることはちゃんと守られている。

この通りベイズ化は期待通りの効果を発揮してくれた様子だが、別の問題も生じている。
N=100, \alpha=1.0 における推定値の平均 \mathbb{E}[\hat\pi_i] を見てみよう。

\pi_i 0.1 0.2 0.3 0.4
\mathbb{E}[\hat\pi_i] 0.18 0.22 0.27 0.32
\rm{std}[\hat\pi_i] 0.09 0.11 0.12 0.13
95%区間 0.06~0.41 0.07~0.48 0.09~0.54 0.11~0.59

分散が \pi_i ごとに異なり、また最尤推定解のそれより小さいなど嬉しいところもあるが、今見てほしいのは最尤推定では成立していた \mathbb{E}[\hat\pi_i]=\pi_i| が成立しなくなっている点である。
これはディリクレ事前分布によって正則化の効果が生じて、解が平準化されているためだ(0.1,0.2 の推定解は大きく、0.3,0.4 のは小さくなっている)。
N=10000 ではまたちゃんと \mathbb{E}[\hat\pi_i]=\pi_i| が成立しており、データが少ないと事前分布に引っ張られるという正しくベイズらしい振る舞いになっている。

では事前分布による正則化の効果を下げるために \alpha を小さくしてみたらどうだろう。N=100,1000,10000 それぞれについて事前分布のパラメータを \alpha=0.1 にしてみたときの推定解の分布が以下の図だ。

f:id:n_shuyo:20210120155522p:plain:w240f:id:n_shuyo:20210120155532p:plain:w240f:id:n_shuyo:20210120155542p:plain:w240

N=10000 では事前分布の違いを感じさせないが、N=100,1000 では \pi_1 が高い確率で 0 に縮退することがわかる。
N=100, \alpha=0.1 における推定値の平均や分散も確認する。

\pi_i 0.1 0.2 0.3 0.4
\mathbb{E}[\hat\pi_i] 0.13 0.19 0.28 0.4
\rm{std}[\hat\pi_i] 0.18 0.23 0.27 0.3
95%区間 0~0.69 0~0.81 0~0.89 0~0.93

\alpha を小さくしたことで、たしかに平準化の効果は押さえられ、\mathbb{E}[\hat\pi_i] は真の割合に近づいている。しかし 95% 区間を見れば分かる通り、一番大きい \pi_4 でさえも 0 に縮退する可能性が十分ある( PRML に書いてあった「関連度自動決定(ARD)」と同じ効果)。
0 に潰れる可能性があるということは、広い範囲の値を取りうるということで分散も大きくなっている。

このようにベイズ推論を用いる場合は事前分布をうまく決める必要がある、という至極当たり前な結論がでたところで*3、Randomized Response のグラフィカルモデルをもう一度よく見てみると、LDA の簡易版(ドキュメントが1個で、単語の生起確率が given)であることがわかる。
となるとせっかくだから LDA ではポピュラーな Collapsed Variational Bayes 推論も試してみたいよね! ということで次回記事に続く。

shuyo.hatenablog.com

Reference

  • Oh, Man-Suk. "Bayesian analysis of randomized response models: a Gibbs sampling approach." Journal of the Korean Statistical Society 23.2 (1994): 463-482.

*1:本稿では触れないが、local differential privacy における privacy budget が一定の範囲で推定値の分散が一番小さくなるのは(式1)の形で表現できることが示せる。

*2:\boldsymbol{X}_{-n} に依存しないので、観測データのループを書かなくていい。今回の Python 実装では、numpy.random.multinomial を選択肢の種類数である M 回呼ぶだけで済んでおり、データ数の N を増やしても実行時間は変わらない。

*3:今回はまだそこまで手を伸ばしてないが、経験ベイズ的アプローチでエビデンスが最大となるようにハイパーパラメータを決定してもおもしろい結果が出るかも?