Mixtureモデルの周辺尤度をモンテカルロ積分で近似してみる
はじめに
この記事では、Mixtureモデルの周辺尤度をモンテカルロ積分で近似するアルゴリズムChib's methodを試します。 Mixtureモデルの周辺尤度は簡単には求まらない(数値的にも)という印象だったので、サンプリングを使って計算する手法がどこまでうまく行くのか検証しようと思いました。
経緯
Non-local prior というモデル選択のための事前分布の紹介qiitaを読んで、興味を持ちNon-local priorの論文も読みました。この論文はクラスタのダブりや空のクラスタが現れる事前確率を0にすることで、周辺尤度を使ったモデル選択(この場合はクラスタ数選択)の性能を上げるというものでした。この論文は面白かったのですが、周辺尤度の計算には疑問が残りました。見てみると周辺尤度の計算にはこの論文(Marin & Robert, 2008)の方法を使っており、この方法はこの論文内でChib's methodと呼ばれる方法をベースとしています。そこで、Chib's methodからMarin & Robertの方法からNon-local priorまですべて計算機実験でためそう!となったわけです。
そして、Chib's methodの検証までで飽きたのでそこまでを記事にすることにしました。
想定するモデル
想定するモデルは以下のようなものです。
- Mixtureモデル
- サンプルはi.i.d.
- パラメタの事前分布はconjugate
Chib's method
まず、Chib's methodを説明します。ベイズの定理より、
が成り立ちます。ここで、はデータ、はパラメタを指します。今回欲しいのは周辺尤度なので、事後分布さえ計算できればOKです。しかも左辺はに依存しないので、として任意の値を代入した確率密度が計算できればよいです。簡単簡単。と、思いきや、肝心なことを忘れています。そうですクラスタを示す潜在変数です。に含めて考えることもできますが、潜在変数とパラメタの同時分布の確率密度を得る必要が出てきて苦しいです。
そこで、潜在変数を周辺化して上式を得ることを考えます。とは簡単です。前者はの事前分布の割合で按分するだけで、後者はそもそも必要ないです。問題はです。簡単には求まりません。
解決法は単純です。MCMCによってのサンプルを得てそれを用いて近似して評価します。すなわち、
と近似します。ここで、はMCMCによる番目のの事後サンプルを示し、はサンプリング回数を示します。 これによって周辺尤度が近似できます。この方法がChib's methodです。 また、Marin & Robertの方法ではMCMCにおけるlabel switchingの問題に対処するためを各クラスタについてパーミュテーションした値全通りで計算して平均を取っています。 今回はとして各クラスタで同じ値を使っています。いい値を選ばないと精度が悪くなるとのことですが、今回は気にしていません。
検証実験
ここからは検証実験の内容になります。
モデル
用いるモデルは超簡単なGaussian mixture modelです。具体的な生成過程は、
です(がになってます)。厳密には記号の定義をするべきですが、しません。ハートで理解してください。
正解のための解析的な計算
ここからが問題です。Chib's methodを検証するためには正確な周辺尤度を計算しないといけません。しかし今回対象にできるモデルはmixtureモデルのみなのでそれができたら苦労しない...。なので、サンプルサイズを小さくし、による積分を厳密に計算できるデータで検証を行いました。が計算できるなら正確な周辺尤度を出せるはず...!つまりは、上の周辺事後分布の
の部分を計算します。 まずです。こちらは、単峰分布の場合の事後分布の計算とほぼ変わりません。すなわち、
のように計算できます。右辺はクラスタkに属するデータについてだけ(クラスタの分け方すなわちは与えられているので可能)単峰の場合の事後分布を計算することになるので、が共役事前分布であることを使って正規化し左辺を求められます。 今回の場合では、
となります。ただし、はクラスタkに属するデータ数です。 を計算するコードは以下になります。
def calc_ln_alpha(x, z, mu_star, D, K): ln_alpha = 0 for k in range(K): x_k = x[z==k] N_k = x_k.shape[0] if N_k == 0: ln_alpha += multivariate_normal.logpdf( mu_star, mean=np.zeros(D), cov=np.identity(D)) continue post_mu_ex = x_k.sum(0) / (N_k + 1) post_mu_cov = np.identity(D) / (N_k + 1) ln_alpha += multivariate_normal.logpdf( mu_star, mean=post_mu_ex, cov=post_mu_cov) return ln_alpha
次にです。ベイズの定理から、
が成り立ちます。ここで、
となり、kに所属するxについてだけ注目したときの、周辺尤度を計算すればよいことがわかります。これは単峰の場合の周辺尤度なので計算できます。具体的には、単峰を考えてを省略するとして、
def calc_ln_betas(x, N, D, K): ln_betas = [] for i in range(pow(K, N)): this_z = [] cluster_str = None if K == 1: cluster_str = ''.zfill(N) else: cluster_str = np.base_repr(i, K).zfill(N) for s in cluster_str: this_z.append(int(s)) this_z = np.array(this_z) ln_beta = 0 for k in range(K): x_k = x[this_z==k] N_k = x_k.shape[0] if N_k == 0: continue ln_betak = - D * N_k * np.log(2 * np.pi) / 2 ln_betak -= D * np.log(N_k + 1) / 2 x_term = 0 for i in range(N_k): x_term += x_k[i] @ x_k[i] x_term -= x_k.sum(0) @ x_k.sum(0) / (N_k+1) x_term /= -2 ln_betak += x_term ln_beta += ln_betak ln_betas.append(ln_beta) ln_betas -= logsumexp(ln_betas) return ln_betas最後に、掛け合わせてKN通りの和を取って、を得ます。 コードは以下です。としていくつかの値を計算して同じ値になったので実装に問題はない...と思ってます。
def calc_true_log_posterior(x, mu_star, N, D, K): ln_etazs = [] ln_betas = calc_ln_betas(x, N, D, K) for i in range(pow(K, N)): this_z = [] cluster_str = None if K == 1: cluster_str = ''.zfill(N) else: cluster_str = np.base_repr(i, K).zfill(N) for s in cluster_str: this_z.append(int(s)) this_z = np.array(this_z) ln_etaz = calc_ln_alpha(x, this_z, mu_star, D, K) + ln_betas[i] ln_etazs.append(ln_etaz) return logsumexp(ln_etazs)
Chib's methodのためのギブスサンプリング実装
こちらは普通にギブスサンプリングです。初期値から、を使った尤度と事前確率の値を正規化した分布からをサンプリング、をgivenとして共役事前分布を利用して計算した事後分布からをサンプリングします。まあ定式化はしなくていいでしょう。ギブスサンプリングとそれを使って対数事後確率を近似するコードは以下です。
def GMM_gibbs_sample(x, K, T): N, D = x.shape mu_sample = np.random.normal(0, 1, size=(K, D)) z_samples = [] mu_samples = [] mu_samples.append(mu_sample) for t in tqdm(range(T)): z_sample = np.zeros((N), dtype=int) for i in range(N): ln_posterior = np.zeros(K) for k in range(K): ll = multivariate_normal.logpdf(x[i], mean=mu_sample[k], cov=np.identity(D)) ln_prior = np.log(1 / K) ln_posterior[k] = ll + ln_prior ln_posterior -= logsumexp(ln_posterior) z_sample[i] = np.random.choice(range(K), p=np.exp(ln_posterior)) z_samples.append(z_sample) for k in range(K): x_k = x[z_sample==k] N_k = x_k.shape[0] mu_sample[k] = np.random.multivariate_normal( x_k.sum(0) / (N_k + 1), np.identity(D) / (N_k + 1)) mu_samples.append(mu_sample.copy()) return np.array(z_samples), np.array(mu_samples)def calc_GS_approx_log_posterior(x, mu_star, D, K): T = 1000 z_samples, _ = GMM_gibbs_sample(x, K, T) n_burn_in = 100 ln_alphas = [] for z_sample in z_samples[n_burn_in:]: ln_alphas.append(calc_ln_alpha(x, z_sample, mu_star, D, K)) lpos_approx = logsumexp(ln_alphas) - np.log(T-n_burn_in) return lpos_approx
カーネル密度推定
加えて、共役事前分布を利用しないケースも考えてみます。共役事前分布でないとを解析的に求められないため、MCMCののサンプルからカーネル密度推定で確率密度関数の値を近似することにします。 この場合のMCMCはギブスサンプリングを使ってしまうとズルなので(上記のギブスサンプリングは共役事前分布を利用している)、Stanでサンプリングすることにします。 Stanコード(を格納する文字列変数を定義するpythonコード)は以下です。
stan_code = """ data { int<lower=0> N; int<lower=0> D; int<lower=0> K; vector[D] X[N]; } parameters { vector[D] mu[K]; } transformed parameters{ vector[K] ps[N]; for (i in 1:N){ for (k in 1:K){ ps[i, k] = normal_lpdf(X[i] | mu[k], 1) + log(1.0/K); } } } model { for (k in 1:K){ mu[k] ~ normal(0, 1); } for (i in 1:N){ target += log_sum_exp(ps[i]); } }このStanコードを使って近似したを得るコードは以下です。
def calc_KDE_approx_log_posterior(x, mu_star, N, D, K): stan_data = {"N": N, "D": D, "K": K, "X": x, } random_seed = np.random.randint(0, 10e6) posterior = stan.build(stan_code, data=stan_data, random_seed=random_seed) fit = posterior.sample(num_chains=4, num_samples=10000) # クラスタを区別しないmuの事後サンプル mu_nuts = fit["mu"] mu_post_samples = mu_nuts.transpose((2, 0, 1)).reshape(-1, D) kde = KernelDensity(kernel='gaussian', bandwidth='silverman').fit(mu_post_samples) log_density = kde.score(mu_star.reshape(1, -1)) return log_density
対数周辺尤度の計算
これで、三者の方法による事後分布の確率密度関数の値の計算ができました。あとは対数周辺尤度の計算です。 事前分布と尤度は共通です。以下のコードで計算できます(として各クラスタで同じ値を用いる前提のコードです)。
def calc_log_prior(mu_star, D, K): lpri = 0 for k in range(K): lpri += multivariate_normal.logpdf(mu_star, np.zeros(D), np.identity(D)) return lpri対数周辺尤度は、事後分布を計算する関数を引数として、以下のように書きました。def calc_log_likelihood(x, mu_star, D): lls = multivariate_normal.logpdf(x, mu_star, np.identity(D)) return lls.sum()
def calc_lml(x, mu_star, N, D, K, calc_log_posterior): ll = calc_log_likelihood(x, mu_star, D) lpri = calc_log_prior(mu_star, D, K) lpos = calc_log_posterior(x, mu_star, N, D, K) return ll + lpri - lpos
実験結果
では、いよいよ三者の比較を行います。確認です。 True: 正解周辺尤度を厳密に計算。共役事前分布が必要で、Nは30が限度。 GS: Chib's methodを使った近似計算。共役事前分布が必要。 KDE: KDEを使ったChib's method。どんなモデルでも使える。KDEのためのハイパーパラメタは必要。
データは以下のコードで生成しました。まあ正解パラメタを当てたいわけではないのでサンプルサイズと次元さえ決められればどんな作り方でも良いです。
def generate_data_and_param(N, D, K): mu = np.random.multivariate_normal(np.zeros(D), np.identity(D), size=K) z = np.random.choice(range(K), size=N) x = [] for i in range(N): x.append(np.random.multivariate_normal(mu[z[i]], np.identity(D))) x = np.array(x) return x, z, muサンプルサイズは3, 5, 10、サンプルの次元は1, 3, 5、クラスタ数は1, 2の全ての組み合わせで作成したデータで実験を行いました。 比較実験は以下のコードで行いました。各条件で生成したデータ一つに対して、各近似手法で100回ずつ計算を行って近似のブレも調べました。
T = 100 settings = [(N, D, K) for N in [3, 5, 10] for D in [1, 3, 5] for K in [1, 2]] for N, D, K in settings: output_filename = f'N{N}D{D}K{K}.pickle' print(output_filename) if os.path.exists(output_filename): continue x, z, mu = generate_data_and_param(N, D, K) mu_star = x.mean(0) true_lml = calc_lml(x, mu_star, N, D, K, calc_true_log_posterior) GS_lmls = [] KDE_lmls = [] for t in tqdm(range(T)): GS_lml = calc_lml(x, mu_star, N, D, K, calc_GS_approx_log_posterior) KDE_lml = calc_lml(x, mu_star, N, D, K, calc_KDE_approx_log_posterior) GS_lmls.append(GS_lml) KDE_lmls.append(KDE_lml) result = [true_lml, GS_lmls, KDE_lmls] with open(output_filename, mode='wb') as f: pickle.dump(result, f)そして可視化はこちらのコードで行いました。
fig, axes = plt.subplots(3, 6, figsize=(30, 15)) fontsize = 20 settings = [(N, D, K) for N in [3, 5, 10] for D in [1, 3, 5] for K in [1, 2]] for t, (N, D, K) in enumerate(settings): i = t // 6 j = t % 6 ax = axes[i][j] condition = f'N{N}D{D}K{K}' output_filename = f'{condition}.pickle' with open(output_filename, mode='rb') as f: true_lml, GS_lmls, KDE_lmls = pickle.load(f) approx_lml_df = pd.DataFrame([GS_lmls, KDE_lmls], index=['GS', 'KDE']).T sns.violinplot(x='variable', y='value', data=pd.melt(approx_lml_df), scale='width', ax=ax) plt.sca(ax) plt.hlines(true_lml, -0.5, 1.5, color='blue') plt.title(condition, fontsize=fontsize) plt.xticks(fontsize=fontsize) plt.yticks(fontsize=fontsize) plt.xlabel('Method', fontsize=fontsize) plt.ylabel('Log marginal likelihood', fontsize=fontsize) plt.plot() plt.tight_layout()結果は以下になります。
まず、肝心のGS(Chib's method)ではそこそこ当たっていそうな感じです。あと、サンプルサイズや次元を大きくすると急激に精度が悪化するようなこともなさそうです(K=2までしか調べてないですし、各条件でデータセットも一つしか試していないので雰囲気です)。KDEだとバイアスが大きそうですね。やり方をちゃんと吟味すればうまく行くかもしれません。 周辺尤度の近似精度としてはそんな感じです。
やり残したこと
近似精度の解釈
ただ、周辺尤度を計算したいときというのはモデル選択を行いたいときが多いはずで、他のモデルでの周辺尤度との差も考慮して誤差を評価したいところです。例えば、近似値のほとんどが真の値の±0.1くらいに収まったとしても、比較したいモデルの真の値に0.001くらいの差しかなければ近似手法の有用性は低そうです(ノイズに差が埋もれてしまう)。
渡辺ベイズ理論における周辺尤度の近似
渡辺ベイズ(ベイズ統計の理論と方法、5.1.4 自由エネルギーの近似)では周辺尤度は逆温度を少しずつ変えた複数回のMCMCにより計算できる一方で、ナイーブなモンテカルロ積分ではうまくいかないと書かれています。 WBICも含めて、渡辺ベイズ理論の手法とも比較したいところでした。畜生、根気さえあれば...!
の吟味
元論文では最尤推定などで事後分布の大きいところを占めるの値をとして設定することが推奨されていました。今回は各クラスタで同じ値を使ってしまったので、理論の上では同じことだとしても、性能に差が出る可能性があります。
高次元データでの性能
今回は低次元データのみでの検証でした。高次元データではどうなるのかわかりません...ので、軽くやってみたところ... ん...?あれ...?若干分布が怪しいですね。中央値を見ればそこそこ近そうですが... さらなる調査が望まれる。
まとめ
いかがでしたか? タイトルの内容を試して、理論値(と言っていいのか?正解の値?)と比較しました。 やり残したことはありますが、個人的には満足です。 では。