Variational Inferenceでガウス分布のパラメータ推定
1変数ガウス分布の推定式の導出
- http://statmodeling.hatenablog.com/entry/variational-bayesian-inference-1
- http://statmodeling.hatenablog.com/entry/variational-bayesian-inference-with-sympy
上の記事を読んだ。variational inferenceという手法があって, これを使うとMCMCサンプリングをせずに事後分布を推定できるのだそうだ. これはとても素晴らしいことだと思う. 参考文献としてPattern Recognizing and Machine Learningの10章を読むといいらしいが高額な書籍で手が出せない. 幸いにして日本語のレビュー記事やスライド資料がたくさんあるので, 別に本がなくても何とかなる.
variational inferenceはパラメータの事後分布を以下で示されるで近似する. 証明は難しいので省略. 上の記事の1. に書いてある.
, はそれぞれ測定データと分布のパラメータである.
今回はvariational inferenceで1変数ガウス分布のパラメータ推定を行う.
まずは何はなくとも尤度関数である.
続いて平均とprecision/精度の事前分布である. は事前分布のパラメータである. precisionは分散の逆数なのだが, これを使うのは共役事前分布がガンマ分布になってくれるからだと思う. 分散を使うと逆ガンマ分布になるので, これが嫌なのだろう.
(1), (2), (3), (4)からを求める. 対数を取ることで掛け算が足し算になるので, に関係ない項を次々にconstにまとめることができる. 共役事前分布を使っているので, 過去の例からしてについてはガウス分布, についてはガンマ分布になるのだろうと予想できる. 手計算による規格化が不要なのでconstは未来永劫に無視できる.
がの二次式なので, はガウス分布 となる. 各パラメータはのの各係数から下のように得られる.
続いてを求める. 以下のようにの1次と対数となる.
これはガンマ分布( )であり, パラメータ は以下となる.
できた.
とりあえず実験
[tex:q^{}(\mu)]と[tex:q^{}(\tau)]のいずれもお互いの期待値を含んでいる. こういう場合は連立方程式として解くか, 適当な初期値を与えてとを順繰りに算出するのを適当な回数こなす方法があるという. 数値的に解く場合は, 収束することは証明されているらしい.
実際に実験してみる. まず準備
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
事後分布のパラメータを更新する関数を用意する.
def update_mu(N, meanX, meanX2, mu0, lambda0, a0, b0, muN, lambdaN, aN, bN): lambdaN = aN / bN * (lambda0 + N) muN = (lambda0 * mu0 + N * meanX) / (lambda0 + N) return lambdaN, muN def update_tau(N, meanX, meanX2, mu0, lambda0, a0, b0, muN, lambdaN, aN, bN): aN = (N + 1) / 2 + a0 bN = b0 + (N + lambda0) / 2 * muN * muN \ - (N * meanX + lambda0 * mu0) * muN \ + (N + lambda0) / (2 * lambdaN) + (N * meanX2 + lambda0 * mu0 * mu0) / 2 return aN, float(bN)
とりあえずStanModelingと同じ条件でを推定してみる. の初期値を3, の初期値を10として, パラメータを10回更新する.
X = pd.DataFrame([1.1, 1.0, 1.3]) mu0 = 0 lambda0 = 1 a0 = 1 b0 = 1 N = X.size meanX = float(X.mean()) meanX2 = float((X*X).mean()) muN = 3 lambdaN = 10 aN = None bN = None for i in range(10): aN, bN = update_tau(N, meanX, meanX2, mu0, lambda0, a0, b0, muN, lambdaN, aN, bN) lambdaN, muN = update_mu(N, meanX, meanX2, mu0, lambda0, a0, b0, muN, lambdaN, aN, bN) print(muN, lambdaN, aN, bN)
0.8500000000000001 1.0958904109589043 3.0 10.95
0.8500000000000001 3.603603603603604 3.0 3.3299999999999996
0.8500000000000001 5.825242718446603 3.0 2.0599999999999996
0.8500000000000001 6.492335437330929 3.0 1.8483333333333332
0.8500000000000001 6.618660946836219 3.0 1.8130555555555554
0.8500000000000001 6.640194697066736 3.0 1.8071759259259257
0.8500000000000001 6.643797285578193 3.0 1.8061959876543208
0.8500000000000001 6.644398097084033 3.0 1.8060326646090534
0.8500000000000001 6.644498242899829 3.0 1.8060054441015088
0.8500000000000001 6.644514934162626 3.0 1.8060009073502514
収束したらしい. 可視化してみるとStanModelingと同じ結果が得られたらしいことが確認できた.
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) q_mu = lambda x: sp.stats.norm.pdf(x, muN, np.sqrt(1/lambdaN)) q_tau = lambda x: sp.stats.gamma.pdf(x, a=aN, scale=1/bN) f = lambda x,y: q_mu(x) * q_tau(y) Z = f(X, Y) cs = plt.contour(X, Y, Z) plt.clabel(cs, inline=1, fontsize=10) plt.xlabel(r'$\mu$') plt.ylabel(r'$\tau$') plt.savefig("vi1_01.png")
データ数を増やして実験
適当な1変数正規分布から50点サンプリングしてを用意する. ここから分布のパラメータを推定し, 事後予測分布を作成してのヒストグラム, 事前分布と比較した. 追記:予測分布はサンプリングした.
# data np.random.seed(0) X = pd.DataFrame(np.random.normal(1, 1, 50), columns=['sample']) N = X.size meanX = float(X.mean()) meanX2 = float((X*X).mean()) # prior distribution mu0 = 0 lambda0 = 1 a0 = 1 b0 = 1 # initial values muN = 3 lambdaN = 10 aN = None bN = None # estimate parameters for i in range(10): aN, bN = update_tau(N, meanX, meanX2, mu0, lambda0, a0, b0, muN, lambdaN, aN, bN) lambdaN, muN = update_mu(N, meanX, meanX2, mu0, lambda0, a0, b0, muN, lambdaN, aN, bN) # prior distribution mu = np.random.normal(mu0, np.sqrt(1/lambda0), 2000) tau = np.random.gamma(a0, 1/b0, 2000) priorX = pd.DataFrame(np.random.normal(mu, np.sqrt(1/tau)), columns=['prior dist.']) # estimate posterior predicitive distribution mu = np.random.normal(muN, np.sqrt(1/lambdaN), 2000) tau = np.random.gamma(aN, 1/bN, 2000) postX = pd.DataFrame(np.random.normal(mu, np.sqrt(1/tau)), columns=['posterior\npredictive dist.']) # plot ax = X.plot.hist(normed=True, bins=15) priorX.plot.kde(ax=ax) postX.plot.kde(ax=ax) ax.set_xlim(-10, 10) ax.set_xlabel('x') ax.set_ylabel(r'Density') plt.savefig("vi1_02.png");
それらしいように見える.
1変数の1ガウス分布の推定ができたので, 次は多変数か混合分布についてやってみたい.
追記: 予測分布が解けない
上では予測分布をサンプリングしたが試しに解けるかやってみた.
予測分布は次の積分となる. にわかには解けない.
とりあえずを解く.
を戻す.
ここで詰まった. これ以上はなぜか解けない. といって数値計算で積分を解く方法は知らない. これ以上は無理.
おわり.