1変数ガウス分布(尤度関数)とNormal-gamma分布(事前分布)から事後分布と事後予測分布を求める
1変数ガウス分布(尤度関数)とNormal-gamma分布(事前分布)から事後分布をvariational inferenceで解くという行為について
前回の1変数正規分布モデルのvariational inferenceでは事後予測分布が導出できなかった. 途中で計算を間違えたのかもしれないし, 解析的に解けない問題だったのかもしれない.
ところであの問題だが、1変数ガウス分布, Normal-gamma分布をパラメータの事前分布とした場合のパラメータの事後分布や事後予測分布は解析的に解けるのだそうである. 解けるところを試しにvariational inferenceしてみて結果を比較したりする, というのが例題の作為であったようだ. そういう基礎的な背景があったようだ. 知らなかった.
というわけで, 普通に解いてみる.
事後分布
こちらが尤度と事前分布になります.
との事後分布を解く. 例によって対数にして考える.
解けた.
ちなみに上のような, の事前分布/事後分布の形をNormal-gamma分布という.
事後予測分布を解く準備
続いて事後予測分布を解くのだが, その前にMurphy, K. P. (2007). "Conjugate Bayesian analysis of the Gaussian distribution."を参考にを解く. これを解いておくと事後予測分布を解くのが楽になる. また, ここでのの解き方は非常に参考になる.
まず, Normal-gamma分布とガウス分布について, 正規化項を除いたとに依存する部分を表す関数とを導入する. また, Normal-gamma分布の正規化項としてを導入する. ガウス分布の正規化項はである.
以下のようにを導出する. このとき, やに依存する項は割って1にならなければならないところがミソである(3行目).
ベイズ統計で尤度と事前分布の積の対数をとって, 定数項=正規化係数を無視して式形からパラメータの事後分布を求めるのは定番のテクニックだと思う. 初めて見たときは感動すら覚えた.
これを逆に考えると, 正規化係数のみに着目すればベイズの定理の右辺の分母 が簡単に導出できるのである. たしかに言われてみればその通りであるがすごい.
事後予測分布
それではある一つのデータについての事後予測分布を解く.
が出てきた. これらはが簡単に示されるのでここから次のように得られる.
つづき.
ただし,
が解けた. これがなにかというと, 自由度, 精度, 平均のStudentのt分布である.
事後予測分布も解けた. すごい.
何がすごいかというと, 面倒くさそうな積分 を解くことなく事後予測分布が得られたことである.
パラメータの事後分布の可視化
可視化してみる, 前回のVariational inferenceの結果とも比較する.
import numpy as np import scipy as sp from scipy import stats import pandas as pd pd.set_option('display.width', 200) import matplotlib.pyplot as plt plt.style.use('ggplot') import base64 from itertools import chain import tempfile import os import subprocess import warnings warnings.filterwarnings('ignore') import sympy import sympy.stats
# test data X = pd.DataFrame([1.1, 1.0, 1.3]) N = X.size meanX = X.mean()[0] # parameters of prior mu0 = 0 lambda0 = 1 a0 = 1 b0 = 1 # parameters of posterior muN = (N * meanX + lambda0 * mu0) / (lambda0 + N) lambdaN = lambda0 + N aN = N / 2. + a0 bN = b0 + 1./2. * X.apply(lambda x: np.power(x - meanX, 2.)).sum()[0] \ + N * lambda0 * np.power(meanX - mu0, 2.) / 2. / (lambda0 + N) print(muN, lambdaN, aN, bN)
0.8500000000000001 4 2.5 1.5050000000000003
# grid delta = 0.05 x = np.arange(-1.0, 3.0, delta) y = np.arange(0.0, 6.0, delta) X, Y = np.meshgrid(x, y) # posterior distribution q_mu = lambda x, y: sp.stats.norm.pdf(x, muN, np.sqrt(1./(lambdaN*y))) q_tau = lambda y: sp.stats.gamma.pdf(y, a=aN, scale=1/bN) f = lambda x,y: q_mu(x, y) * q_tau(y) Z = f(X, Y) # posterior distribution with Variational inference muN_vi, lambdaN_vi, aN_vi, bN_vi = 0.8500000000000001, 6.644514934162626, 3.0, 1.8060009073502514 q_mu_vi = lambda x: sp.stats.norm.pdf(x, muN_vi, np.sqrt(1/(lambdaN_vi))) q_tau_vi = lambda x: sp.stats.gamma.pdf(x, a=aN_vi, scale=1/bN_vi) f_vi = lambda x,y: q_mu_vi(x) * q_tau_vi(y) Z_vi = f_vi(X, Y) # plot cs = plt.contour(X, Y, Z) cs_vi = plt.contour(X, Y, Z_vi, colors='red', alpha=0.5) plt.clabel(cs, inline=1, fontsize=10) plt.xlabel(r'$\mu$') plt.ylabel(r'$\tau$') plt.title('posterior distribution') lines = [cs.collections[0], cs_vi.collections[0]] labels = ['Conjugate Bayesian analysis','Variational inference'] plt.legend(lines, labels) plt.savefig("1.png")
事後予測分布の可視化
解析的に得られた事後予測分布とサンプリングによる事後予測分布を比較する.
precision = aN * lambdaN / bN / (lambdaN + 1.) # grid delta = 0.1 df = pd.DataFrame() df['X'] = np.arange(-4.0, 6.0, delta) # posterior predicitive distribution df['calculated'] = sp.stats.t.pdf(df.X, df=2.*aN, loc=muN, scale=np.sqrt(1/precision)) # sampled posterior predictive distribution tau = np.random.gamma(aN, 1/bN, 2000) mu = np.random.normal(muN, np.sqrt(1./lambdaN/tau)) sampled = pd.DataFrame(np.random.normal(mu, np.sqrt(1/tau)), columns=['sampled']) # plot ax = df.plot(x='X', y='calculated') sampled.plot.kde(ax=ax) ax.set_xlim(-4, 6) ax.set_xlabel('x') ax.set_ylabel(r'Density') ax.set_title('posterior predictive distribution') plt.savefig("2.png");
だいたい一致. よさそう.