多変量正規分布+Normal-Wishart分布(事前分布)の事後分布
前回1変数をやったので多変量もやりたい.
導出
分布の定義
まずは多変量正規分布である. は の次元数.
つづいてWishart分布である. について
統計モデル
データ は多変量正規分布に従う. 事前分布をNormal-Wishart分布とする.
が事前分布のパラメータである.
との事後分布
例によって対数をとって考える.
ひとまずの項であるをやっつけるとする.
をやっつける.
を対数事後分布に戻し, Normal-Wishart分布になるよう変形.
事後分布が導出できた. そんなに難しくない.
の周辺事後分布
との事後分布が得られたわけだが, ここからの周辺分布は導出できるのか考えてみる. ちなみにの周辺事後分布は明らかに事後分布のWishart分布の部分である.
要するに下のNormal-Wishart分布のの平均に対する周辺分布を導出して, パラメータを事後分布のパラメータに置き換えればよい.
Wishart分布の正規化項を導入すると計算しやすい.
しこしこ解く.
について積分し, の周辺分布を求める.
ここでの処理は光成滋生 (2017). "パターン認識と機械学習の学習 普及版"のp79が参考になる.
よって以下のように導出できる.
これは, 精度, 平均 , 自由度のStudent's t分布である.
そういうわけで, は, 精度, 平均 , 自由度のStudent's t分布となる.
解けた.
事後予測分布
任意のについての事後予測分布 を導出する.
まず を導出する. 前回やったように, 正規化項以外の割り算は割って1にならなければならない.
は以下となる.
による事後分布のパラメータが出てきた. このパラメータを添字 で表現すると以下となる.
は, 精度 , 平均 , 自由度 のstudentのt分布となる.
のStudent's t分布の精度を で割った分布なわけだがそこに意味があるのかは分からない.
なんにせよ導出は全部できた.
可視化と検証
計算して可視化してみる. Rstanでサンプリングもしてみて結果を比較してみる. 一致しなければ上の導出が誤っていることになる.
import numpy as np import scipy as sp from scipy import stats from scipy import special import pandas as pd pd.set_option('display.width', 200) import matplotlib.pyplot as plt import matplotlib plt.style.use('ggplot') import os import subprocess import warnings warnings.filterwarnings('ignore')
テストデータと事前分布
適当な二次元データをサンプリングする.
# 適当な二次元データ np.random.seed(0) mu = np.array([5, 3]) cov = np.array([[5,2],[2,1]]) N = 100 X = pd.DataFrame(sp.random.multivariate_normal(mu, cov, N), columns=['x', 'y']) X.plot(kind='scatter', x='x', y='y').set_aspect('equal') plt.savefig('1.png', bbox_inches='tight', pad_inches=0)
事前分布を適当に決める.
# 適当な事前分布のパラメータ mu0 = np.array([0, 0]).reshape(2, 1) # 縦ベクトル lambda0 = 1 W0 = np.array([[0.01, 0], [0, 0.01]]) nu0 = 1.1
事後分布の計算
meanX = np.array(X.mean()).reshape(2, 1) muN = (N * meanX + lambda0 * mu0) / (N + lambda0) lambdaN = N + lambda0 WNinv = np.linalg.inv(W0) \ + (N * lambda0) / (N + lambda0) * (mu0 - meanX) @ (mu0 - meanX).T \ + X.apply( lambda x: (np.array(x).reshape(2, 1) - meanX) @ (np.array(x).reshape(2, 1) - meanX).T , axis=1 ).sum() WN = np.linalg.inv(WNinv) nuN = N + nu0 print(f"muN : {muN.tolist()}") print(f"lambdaN : {lambdaN}") print(f"WN : {WN.tolist()}") print(f"nuN : {nuN}")
muN : [[4.9302018085299], [3.0252711023662453]] lambdaN : 101 WN : [[0.002450141567402565, -0.002569876767106904], [-0.002569876767106904, 0.007340956277110485]] nuN : 101.1
計算はできたがこれが正しいのかはわからない.
RStanで事後分布サンプリング
別の手法で事後分布を求めて比較する必要が有る. というわけでRStanで同じモデルでMCMCサンプリングする. 保存するのはとのサンプリング結果と, 各イテレーションにおいてとをパラメータとする2変量正規分布から1点サンプリングしたものである.
X.to_csv("x.csv", index=False) from subprocess import Popen, PIPE R = Popen(["R", "--vanilla", "--silent"], stdin=PIPE, stdout=PIPE) R.stdin.write(b''' options(width=200, echo=F, digits=4) sink(file="/dev/null") suppressMessages({ library(rstan) }) x <- read.table("x.csv", sep=",", header=T) data = list( N = nrow(x), x = x, mu0 = c(0, 0), lambda0 = 1, W0 = matrix(c(0.01, 0, 0, 0.01), nrow=2, ncol=2), nu0 = 1.1 ) model = " data { int<lower = 0> N; vector[2] x[N]; vector[2] mu0; real lambda0; cov_matrix[2] W0; real nu0; } transformed data { cov_matrix[2] W0inv; W0inv = inverse(W0); } parameters { vector[2] mu; cov_matrix[2] Sigma; } model { target += inv_wishart_lpdf(Sigma | nu0, W0inv); target += multi_normal_lpdf(mu | mu0, Sigma / lambda0); for (n in 1:N) { target += multi_normal_lpdf(x[n] | mu, Sigma); } } generated quantities { vector[2] sampled_x; sampled_x = multi_normal_rng(mu, Sigma); } " fit <- stan( model_code = model, data = data, iter = 20000, warmup = 10000, thin = 1, cores = 4, seed = 1 ) sink() write.table(rstan:::extract(fit), file="stan.csv", row.names=F, sep=",") fit ''') print(R.communicate()[0].decode('ascii')) R.kill()
>
> options(width=200, echo=F, digits=4)
Inference for Stan model: 28d63f23d0d81288800b4248b89bfb14.
4 chains, each with iter=20000; warmup=10000; thin=1;
post-warmup draws per chain=10000, total post-warmup draws=40000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[1] 4.93 0.00 0.26 4.43 4.76 4.93 5.10 5.44 35071 1
mu[2] 3.03 0.00 0.15 2.74 2.93 3.03 3.13 3.32 35489 1
Sigma[1,1] 6.57 0.01 0.95 4.97 5.90 6.48 7.14 8.67 35641 1
Sigma[1,2] 2.30 0.00 0.45 1.53 1.99 2.26 2.58 3.30 32605 1
Sigma[2,1] 2.30 0.00 0.45 1.53 1.99 2.26 2.58 3.30 32605 1
Sigma[2,2] 2.20 0.00 0.32 1.66 1.97 2.16 2.39 2.90 35722 1
sampled_x[1] 4.93 0.01 2.58 -0.12 3.20 4.93 6.66 10.01 39676 1
sampled_x[2] 3.02 0.01 1.49 0.11 2.04 3.02 4.03 5.94 39455 1
lp__ -397.57 0.01 1.60 -401.50 -398.41 -397.24 -396.38 -395.43 17924 1
Samples were drawn using NUTS(diag_e) at Fri Oct 13 01:54:59 2318.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
サンプリング結果のサマリを表示してみると. Rhatがどれも1でn_effも大きいので十分収束したんだと思う. サンプリング結果はCSV経由でPython側で取得する.
stan = pd.read_csv('stan.csv') stan.columns = [f"{c}.stan" for c in stan.columns] stan.columns
Index(['mu.1.stan', 'mu.2.stan', 'Sigma.1.1.stan', 'Sigma.2.1.stan', 'Sigma.1.2.stan', 'Sigma.2.2.stan', 'sampled_x.1.stan', 'sampled_x.2.stan', 'lp__.stan'], dtype='object')
を比較
2次元なので比較が難しい. まずはヒートマップを書いてみる. 多変量のStudent's t分布のPDFは自作する. 多分合ってると思う.
def student_t_pdf(x, m, Lambda, df): D = len(m) m = m.reshape(D, 1) x = x.T return np.diag( sp.special.gamma(D/2. + df/2.) / sp.special.gamma(df/2.) * np.power(np.linalg.det(Lambda), 1/2.) / np.power(np.pi * df, D/2.) * np.power(1 + (x-m).T@Lambda@(x-m)/df, -(D+df)/2.) )
fig, axes = plt.subplots(ncols=3, figsize=(16, 4)) delta = 0.05 gx = np.arange(3.5, 6.5, delta) gy = np.arange(1.5, 4.5, delta) gxx, gyy = np.meshgrid(gx, gy) gxgy = np.c_[gxx.ravel(), gyy.ravel()] zz = student_t_pdf( gxgy, muN, lambdaN * (nuN + 1 - 2) * WN, (nuN + 1 - 2) ) im = axes[0].imshow( zz.reshape(len(gx), len(gy)), interpolation='none', origin='lower', extent=[3.5, 6.5, 1.5, 4.5], cmap=matplotlib.cm.rainbow ) axes[0].set_title('derivated $\mu$ posterior dist.') axes[0].set_xlabel('x') axes[0].set_ylabel('y') fig.colorbar(im, ax=axes[0]) x = stan['mu.1.stan'] y = stan['mu.2.stan'] H = np.histogram2d(x, y, bins=[gx, gy], normed=True) im = axes[1].imshow( H[0].T, interpolation='none', origin='lower', extent=[3.5, 6.5, 1.5, 4.5], cmap=matplotlib.cm.rainbow ) axes[1].set_title('Stan sampled $\mu$ posterior dist.') axes[1].set_xlabel('x') axes[1].set_ylabel('y') fig.colorbar(im, ax=axes[1]) stan.head(1000).plot(kind='scatter', x='mu.1.stan', y='mu.2.stan', alpha=0.1, ax=axes[2]).set_aspect('equal') axes[2].set_xlim(3.5, 6.5) axes[2].set_ylim(1.5, 4.5) axes[2].set_title('Stan sampled posterior $\mu$') axes[2].set_xlabel('x') axes[2].set_ylabel('y') plt.savefig('2.png', bbox_inches='tight', pad_inches=0)
左が の事後分布, 中央がStanがサンプリングした の事後分布, 右はStanのサンプリング結果をそのままプロットしたもの. ヒートマップはだいたい同じ形だけれども, 赤い領域の大きさちょっと違う. Stanの方が若干ばらつきが大きいように見える.
単にStanのサンプルサイズの小ささゆえに2Dヒストグラムがガタついているのが目視での差異の原因かもしれない. だいたい同じなのでこれはこれでよしとしておく.
続いて各次元の周辺事後分布を比較してみると, こちらはぴったり一致している.
delta = 0.05 gridx = np.arange(1, 8, delta) D = 2 data = pd.DataFrame(gridx, columns=['x']) data['mu1'] = pd.DataFrame(sp.stats.t.pdf( gridx, df=(nuN + 1 - D), loc=muN.reshape(2)[0], scale=np.sqrt(np.linalg.inv(lambdaN * (nuN + 1 - D) * WN)[0,0]) )) data['mu2'] = pd.DataFrame(sp.stats.t.pdf( gridx, df=(nuN + 1 - D), loc=muN.reshape(2)[1], scale=np.sqrt(np.linalg.inv(lambdaN * (nuN + 1 - D) * WN)[1,1]) ))
fig, axes = plt.subplots(ncols=2, figsize=(12,4)) data.plot(x='x', y='mu1', ax=axes[0]) stan[['mu.1.stan']].plot.kde(ax=axes[0]) data.plot(x='x', y='mu2', ax=axes[1]) stan[['mu.2.stan']].plot.kde(ax=axes[1]) axes[0].set_xlim(2.5, 7.5) axes[1].set_xlim(1.5, 4.5) plt.savefig('3.png', bbox_inches='tight', pad_inches=0)
まあいいんじゃないですかね.
を比較
これこそどうやってStanのサンプリング結果と比較したらよいのか分からない. だれかおしえてくれ.
簡単に, の事後分布から (分散共分散行列)をサンプリングして各要素の周辺分布をStanのサンプリング結果のSigmaと比較しよう.
sampled_Sigma = pd.DataFrame( sp.stats.invwishart.rvs(df=nuN, scale=np.linalg.inv(WN), size=10000).reshape(10000, 4), columns=['Sigma.1.1', 'Sigma.1.2', 'Sigma.2.1', 'Sigma.2.2'] ) fig, axes = plt.subplots(ncols=3, figsize=(14,4)) sampled_Sigma[['Sigma.1.1']].plot.kde(ax=axes[0]) sampled_Sigma[['Sigma.1.2']].plot.kde(ax=axes[1]) sampled_Sigma[['Sigma.2.2']].plot.kde(ax=axes[2]) stan[['Sigma.1.1.stan']].plot.kde(ax=axes[0]) stan[['Sigma.1.2.stan']].plot.kde(ax=axes[1]) stan[['Sigma.2.2.stan']].plot.kde(ax=axes[2]) plt.savefig('4.png', bbox_inches='tight', pad_inches=0)
一緒ですね. これ以上は深くつっこまない.
解析的に解けるのにStanのサンプリング結果と比較するための事後分布をサンプリングするというのはいかにも筋がわるい.
条件付き事後予測分布の比較
あとは事後予測分布を比較したい. Stanではsampled_xとしてサンプリングしているわけだが.
のときはヒートマップと各次元の周辺事後分布を比較したが, 目的ベースでいくと, の第2成分がある のときの第1成分の事後予測分布の条件付き確率分布を算出してStanと比較したい. これができると様々な条件に対する予測が簡単にできることになる.
そのためには事後予測分布である多変量Student's t分布の条件付き確率を導出しなければならない.
Student's t分布の確率密度関数 は以下である.
ここで
であるときに, が欲しいわけである.
多変量正規分布の条件付き確率によると,
であるとき,
である.
また, Student'sのt分布は以下のように変形できる(PRML 式2.161).
上記より, 以下が得られる.
(多変量正規分布の条件付き確率やPRML 式2.161をリファレンスからそのまま引いているので)簡単じゃないか.
すでに導出したように, 事後予測分布に当てはめると, 上の式のパラメータは, 精度 , 平均 , 自由度 である.
Student's t分布の条件付き確率がわかったので, x1|x2=3の場合の条件付き確率を可視化し, stanのサンプリング結果と比較する. Stanの方は, 2.9 <=x2<=3.1の範囲にあるサンプリング結果だけを抽出してヒストグラムと密度分布を可視化してみる.
# x2 = yのときのxの条件付き確率を算出する def conditional_student_t_pdf(x, y, m, Lambda, df): D = len(m) p = len(x) Sigma = np.linalg.inv(Lambda) Sigma11 = Sigma[0:p, 0:p] Sigma12 = Sigma[0:p, p:] Sigma21 = Sigma[p:, 0:p] Sigma22 = Sigma[p:, p:] mu1 = m[:p, 0] mu2 = m[p:, 0] cmu = mu1 + Sigma12@(Sigma22.T)*(y-mu2) cSigma = Sigma11 - Sigma12@np.linalg.inv(Sigma22)@Sigma21 cLambda = np.linalg.inv(cSigma) return student_t_pdf(x.T, cmu, cLambda, df)
delta = 0.05 gridx = np.arange(-5, 15, delta) D = 2 data = pd.DataFrame(gridx, columns=['x']) # x2=3の場合の各xグリッド点の条件付き確率 data['derivated'] = pd.DataFrame(conditional_student_t_pdf( np.array([gridx]), np.array([[3]]), muN, WN * lambdaN * (nuN + 1 - 2) / (1 + lambdaN), (nuN + 1 - 2) ))
ax = data.plot(x='x', y='derivated') # sampling結果をx2の値の範囲で絞り込み cstan = stan[ (2.9 <= stan['sampled_x.2.stan']) & (stan['sampled_x.2.stan'] <= 3.1) ] # 十分な個数があることを確認, 1次元だと2000は必要と緑本に書いてあった print(cstan.count()) cstan[['sampled_x.1.stan']].plot.hist(normed=True, bins=50, ax=ax) cstan[['sampled_x.1.stan']].plot.density(ax = ax) ax.set_xlim(-5, 15) ax.set_title('conditional posterior predictive dist.(x1|x2=3)') plt.savefig('5.png', bbox_inches='tight', pad_inches=0)
mu.1.stan 2188 mu.2.stan 2188 Sigma.1.1.stan 2188 Sigma.2.1.stan 2188 Sigma.1.2.stan 2188 Sigma.2.2.stan 2188 sampled_x.1.stan 2188 sampled_x.2.stan 2188 lp__.stan 2188 dtype: int64
よく一致しますね. よかったですね. やりましたね.
事後予測分布の比較
最後に事後予測分布のヒートマップを描いて終わる.
fig, axes = plt.subplots(ncols=3, figsize=(16, 4)) delta = 0.5 xmin = -6 xmax = 16 ymin = -6 ymax = 16 gx = np.arange(xmin, xmax, delta) gy = np.arange(ymin, ymax, delta) gxx, gyy = np.meshgrid(gx, gy) gxgy = np.c_[gxx.ravel(), gyy.ravel()] zz = student_t_pdf( gxgy, muN, WN * lambdaN * (nuN + 1 - 2) / (1 + lambdaN), (nuN + 1 - 2) ) im = axes[0].imshow( zz.reshape(len(gx), len(gy)), interpolation='none', origin='lower', extent=[xmin, xmax, ymin, ymax], cmap=matplotlib.cm.rainbow ) axes[0].set_title('derivated posterior predictive dist.') axes[0].set_xlabel('x') axes[0].set_ylabel('y') axes[0].set_yticks(np.arange(-5, 20, 5)) fig.colorbar(im, ax=axes[0]) x = stan['sampled_x.1.stan'] y = stan['sampled_x.2.stan'] H = np.histogram2d(x, y, bins=[gx, gy], normed=True) im = axes[1].imshow( H[0].T, interpolation='none', origin='lower', extent=[xmin, xmax, ymin, ymax], cmap=matplotlib.cm.rainbow ) axes[1].set_title('Stan sampled posterior predictive dist.') axes[1].set_xlabel('x') axes[1].set_ylabel('y') axes[1].set_yticks(np.arange(-5, 20, 5)) fig.colorbar(im, ax=axes[1]) stan.head(10000).plot(kind='scatter', x='sampled_x.1.stan', y='sampled_x.2.stan', alpha=0.1, ax=axes[2]).set_aspect('equal') axes[2].set_xlim(-6, 16) axes[2].set_ylim(-6, 16) axes[2].set_title('Stan sampled posterior x') axes[2].set_xlabel('x') axes[2].set_ylabel('y') axes[2].set_yticks(np.arange(-5, 20, 5)) plt.savefig('6.png', bbox_inches='tight', pad_inches=0)
おしまい.