無限混合正規分布 cont.

前回の補足。

(たぶん)実装可能なメモ書きがこちら

psuedoコード代わりのRコード。
データを入れれば、一応動くけれど、ナイーブな実装なので超遅い(Rでも)。

# 1. prepare dataset
X = ...

# 2. set up hyperparameters
theta0 = list(NULL)
theta0$m = c(0,0)
theta0$xi = 1e-2
theta0$eta = 2*D+2+1e-1
theta0$B = 1e-1*diag(D)*1

# 3. Gibbs sampling
DPM = function(X, theta0){

	# 3.1 functions

	# logarithm of multivariate gamma function
	ln_mgamma = function(x, d){
		0.25 * d * (d - 1) * log(pi) + sum(lgamma(x + 0.5 * (1 - 1:d)))
	}

	# calculate posterior of NIW
	posterior = function(X, theta0){

		if(class(X) == "numeric"){
			X = t(X)
		}

		N = nrow(X)
		D = ncol(X)

		theta = list(NULL)
		theta$xi = theta0$xi + N
		theta$m = (theta0$xi * theta0$m + apply(X, 2, sum)) / theta$xi
		theta$eta = theta0$eta + N
		theta$B = theta0$B + theta0$xi * (theta0$m %*% t(theta0$m)) + t(X) %*% X - theta$xi * (theta$m %*% t(theta$m))

		return(theta)

	}

	# logarithm of the gaussian marginal likelihood for NIW parameter distribution
	ln_marginal_likelihood = function(X, theta0){

		if(class(X) == "numeric"){
			X = t(X)
		}

		N = nrow(X)
		D = ncol(X)

		theta = posterior(X, theta0)

		a = 0.5 * (theta0$eta - D - 1)
		ap = 0.5 * (theta$eta - D - 1)
	
		retval = -0.5*D*N*log(pi) + 0.5*D*log(theta0$xi/theta$xi) + ln_mgamma(ap, D) - ln_mgamma(a, D) + 
				a * determinant(theta0$B)$modulus - ap * determinant(theta$B)$modulus

		return(retval)

	}

	# 3-2. initialize MCMC state & parameters
	N = nrow(X)
	D = ncol(X)
	Z = numeric(N)
	K = 0
	alpha = 1e-1

	itmax = 500

	# 3-3. execute gibbs sampling
	for(it in 1:itmax){

		# randomize samples
		idx = sample(1:N, N)

		# sample z for each x
		for(n in idx){

			Z[n] = 0
			x = t(X[n,])

			lnLn = numeric(K+1)
			if(K > 0){
				for(k in 1:K){
					dat_idx = which(Z == k)
					Nz = length(dat_idx)
					if(Nz > 0){
						lnLn[k] = log(Nz) + ln_marginal_likelihood(x, posterior(X[dat_idx,], theta0))
					}else{
						lnLn[k] = -Inf
					}
				}
			}
			lnLn[K+1] = log(alpha) + ln_marginal_likelihood(x, theta0)
			
			Z[n] = sample(1:(K+1), 1, prob=exp(lnLn - max(lnLn)))

			if(Z[n] == (K+1)){
				K = K + 1
			}

		}

	}	
}

無限混合正規分布

と書くと、とても格好良い。
最終兵器感がある。

別名、正規分布版のDirichlet Process Mixture(DPM)。ノンパラメトリックベイズ推定の1つである。

詳細は、M.D. Escober, M. West, “Bayesian Density Estimation and Inference Using Mixtures”やRadford M. Neal, “Markov Chain Sampling Methods for Dirichlet Process Mixture Models”に詳しい。

連続の確率分布を無限の点で離散化するDirichlet Processを混合モデルの事前分布に置くことで、理論的には無限の混合が可能になる・・・らしい。原理を理解するには、測度論という抽象数学が必要かつ自分にはハードルが高すぎなので、諦めてそういうものかと納得して使っている。

DPMは、実利的には、データの複雑さに応じてモデルの混合数を自動的に推定するための枠組みである。実装の詳細は、別の機会にまとめるが、Neal’s algorithm 3 (collapsed gibbs sampler) を正規分布の尤度関数+Normal-inverse-Wishartのパラメータ分布に対して、実行するとこんな感じになる。

iGM

左下の図が、分割数の学習による変化を表したグラフである。最初は単純な混合モデルに始まり、次第にデータに適した形で複雑さを学習していっていることがわかる。

通常の混合正規分布推定の場合、混合数は予め与えられていなければならない。混合数として、データ数Nを指定しておけば、同じことは一応できるが、全要素について計算しつづけなければならず、負荷は非常に大きい。一方、DPMの場合は、現時点の混合数K+1のみを扱えばよいので、遥かに効率的である。

個人的にはDeep Learningのようにブラックボックスなモデルよりも、こういうわかりやすい(?)モデルのほうが好き。

変分ベイズの定式化から近似事後分布の導出まで(3)

PRML10.1の穴埋め、ラストです。

話の流れは以下になります。

  1. 変分ベイズの定式化
  2. 準備:変分法
  3. 変分事後分布の導出(この記事)

今回は、前々回に定式化した変分ベイズの汎関数最大化問題を、変分法によって解きます。

ここで、たった一つだけ仮定を入れます。それは、近似しようとしている事後分布において、確率変数Z=\{Z_{1}, ..., Z_{L}\}の各Z_{i}が互いに独立であるというものです。手計算できる場合、他の手法により効率的に近似事後分布を発見できる場合、など任意の変数間の非独立性を仮定しても構いません。ここでは、とりあえず全ての変数が独立であるとの仮定から話をしていきます。

また、本来、Fに対して、ラグランジュの未定乗数法により、q(Z)の拘束条件である\int q(Z) dZ = 1を入れるのですが、煩雑になるので、省略します。入れて解いてもただ表記が面倒なだけで、本質的なむずかしさはありません。


3. 変分事後分布の導出

まず、変分ベイズの計算をするにあたって1つの仮定を行う。確率変数Z=\{Z_{1},...,Z_{L}\}は、変分事後分布q(Z)においてそれぞれ独立である。すなわち、

(1)   \begin{equation*} q(Z_{1},Z_{2},...,Z_{L}) = q(Z_{1})q(Z_{2})...q(Z_{L}) \end{equation*}

Fを改めて表記すると

(2)   \begin{eqnarray*} F[q(Z)] &=& F[q(Z_{1}),...,q(Z_{L})] \\ &=& \int q(Z_{1})...q(Z_{L}) \ln \frac{p(D,Z)}{q(Z_{1})...q(Z_{L})} dZ_{1}...dZ_{L} \end{eqnarray*}

ここで、q(Z_{i})以外の変分事後分布を固定して整理すると、(2)はq(Z_{i})についての次のような汎関数となり、前回の変分法で対象としている汎関数の形態に一致する。

(3)   \begin{eqnarray*} F[q_{i}] &=& \int  Q(Z_{i},q_{i})  dZ_{i} \\ Q(Z_{i},q_{i}) &=& \int q_{1}...q_{L} \ln \frac{p(D,Z)}{q_{1}...q_{L}} dZ_{-i} \\ \int q_{i} dZ_{i} &=& 1 \\ q_{k} &\equiv& q(Z_{k}),1 k=1,...,L \\ dZ_{-i} &\equiv& dZ_{1}...dZ_{i-1}dZ_{i+1}...dZ_{L} \end{eqnarray*}

前回の結果から、F[q(Z_{i})]を極大化するq(Z_{i})は次の条件を満たす。

(4)   \begin{eqnarray*} \frac{\partial Q}{\partial q_{i}} = 0 \\ \frac{\partial^{2} Q}{\partial q_{i}^{2}} < 0 \\ \end{eqnarray*}

まず、停留条件から計算する。qを変数とみなして通常の微分を行うと、

(5)   \begin{eqnarray*} \frac{\partial Q}{\partial q_{i}} &=&  \int q_{1}...q_{i-1}q_{i+1}...q_{L} \ln p(D,Z) dZ_{-i}\\ &-& \int q_{1}...q_{i-1}q_{i+1}...q_{L} (\ln q_{1}+...+\ln q_{L}) dZ_{-i} \\ &-& \int q_{1}...q_{i-1}q_{i+1}...q_{L} dZ_{-i} \\ &=& \langle \ln p(D,Z) \rangle_{q_{-i}} - \ln q_{i} + const.\\ &=& 0 \end{eqnarray*}

\langle \ln p(D,Z) \rangle_{q_{-i}}は、\ln p(D,Z)q_{i}を「除く」分布上での期待値である。const.は、最大化しようとしているq_{i}以外の項をまとめている。本来、ここにラグランジュの未定乗数が加わり、正規化項を構成するので任意でよい。以上から、次の関係を得る。

(6)   \begin{equation*} \ln q_{i} = \langle \ln p(D,Z) \rangle_{q_{-i}} + const. \end{equation*}

(7)   \begin{equation*} q_{i} \propto exp\{ \langle \ln p(D,Z) \rangle_{q_{-i}} \} \end{equation*}

(6),(7)は同値で、さらに(6)はPRML10.1の10.9式と同値である。最後に極大条件を確認しておく。

(8)   \begin{eqnarray*} \frac{\partial^{2} Q}{\partial q_{i}^{2}} =  - \frac{1}{q_{i}} < 0\\ \end{eqnarray*}


以上で、当初の目標のPRML10.1の穴埋めである「変分ベイズの定式化から近似事後分布の導出」は完成です。(7)式は独立を仮定した確率変数の数だけ出てくるので、それらを連立して解いていきます。これは完全に問題依存で、紙をたくさん消費してがんばります。

変分ベイズの定式化から近似事後分布の導出まで(2)

前回に引き続きPRML10.1の穴埋めです。

話の流れは以下になります。

  1. 変分ベイズの定式化
  2. 準備:変分法(この記事)
  3. 変分事後分布の導出

前回、

    \[ p(Z|D) \approx \arg \max F[q(Z)] \]

    \[ F[q(Z)] =  \int q(Z) \frac{ \ln p(D, Z) }{ \ln q(Z) } dZ \]

という形で変分ベイズ法による潜在変数・パラメータ諸々込みの確率変数ZのデータDの下での事後分布近似問題の定式化を行いました。

今回は、この最大化問題を解くためのツールとしての変分法の紹介です。そんなに勉強したわけではないので、変な所があれば指摘歓迎です。参考にしたのは、次の2冊です。Calculus of Variationsは必要十分な感じでとても読みやすかったです。日本語だとべんりな変分原理は物理の実例が多めですが、行間もちょっと多めかも?

webリソースでは、こことかここが参考になりました。


2. 準備:変分法

(1)   \begin{equation*} I[y] = \int L(x_{1}, x_{2}, y) dx_{1}dx_{2} \end{equation*}

のような状況を考える。yx_{1},x_{2}の関数であり、Iは関数y(x_{1},x_{2})を引数にとる汎関数と呼ばれる(汎関数の文脈では、yのことを変関数と呼ぶ)。この時I[y]に極値を与える変関数yを求める問題を変分問題といい、通常の関数の極値を求める問題に対応する。ここでは、最大化問題に限定しておく。

VB2

 

ここで、Iに極値を与える変関数をy_{0}と置いて、y_{0}\epsilon \eta(x_{1}, x_{2})だけ微小変化させた関数をyとして再定義する。すなわち、

(2)   \begin{equation*} y(x_{1}, x_{2}) = y_{0}(x_{1}, x_{2}) + \epsilon \eta(x_{1}, x_{2}) \end{equation*}

さらに、yy_{0}の汎関数値の差を\epsilonの関数として次のように定義する。

(3)   \begin{equation*} \Delta I(\epsilon) = I[y] - I[y_{0}] \end{equation*}

(3)式のI[y]\epsilon = 0周り、すなわちyy_{0}と一致する点で、二次までTaylor展開する。

(4)   \begin{equation*} I[y] = I[y_{0}] + \epsilon \frac{dI}{d\epsilon} \biggm|_{\epsilon=0} + \frac{1}{2} \epsilon^{2} \frac{d^{2}I}{d\epsilon^{2}} \biggm|_{\epsilon=0} + O(\epsilon^3) \end{equation*}

\epsilon=0の時、I[y]I[y_{0}]に還元されることに注意する。後は、Taylor展開そのままである。(3)式に(4)を代入すると、

(5)   \begin{eqnarray*} \Delta I(\epsilon) &=& \epsilon \frac{dI}{d\epsilon} \biggm|_{\epsilon=0} + \frac{1}{2} \epsilon^{2} \frac{d^{2}I}{d\epsilon^{2}} \biggm|_{\epsilon=0} + O(\epsilon^3) \\ &=& \delta I + \delta^{2} I + O(\epsilon^3) \end{eqnarray*}

(5)により、y_{0}の極大点付近でのyの微小変化量を表現する関数が構成できた。ここで、\delta Iを第一変分、\delta^{2} Iを第二変分と呼ぶ。

y_{0}においてIが極大であるための条件は次の2つである。

  1. \delta I = 0
  2. \delta I^{2} < 0

\Delta Iは(極大値から少しずれた関数値)から(極大値)を減じているので、\Delta Iは負値にならなければならない。第一変分は\epsilonが正負値を取りうるので、0にならなければならず、条件1が与えられる。これを停留条件と呼ぶ。条件2は第一変分を0とした結果から、第2変分を負にすることが要請される。これがy_{0}が極大であるための条件である。

まず、条件1を詳細化する。

(6)   \begin{eqnarray*} \delta I &=& \epsilon \frac{dI}{d\epsilon} \biggm|_{\epsilon=0} \\ &=& \epsilon \cdot \frac{d}{d \epsilon} \int L(x_{1}, x_{2}, y) dx_{1}dx_{2} \biggm|_{\epsilon=0} \\ &=& \epsilon \cdot \int \frac{d L}{d \epsilon} dx_{1}dx_{2} \biggm|_{\epsilon=0} \\ &=& \epsilon \cdot \int \frac{\partial L}{\partial y} \frac{\partial y}{\partial \epsilon} dx_{1}dx_{2} \biggm|_{\epsilon=0} \\ &=& \int \frac{\partial L}{\partial y} \epsilon \eta dx_{1}dx_{2} \\ &=& 0 \end{eqnarray*}

\epsilon \etaは任意なので、(6)式が恒等的に0になるためには、次の条件を満たせば良い。

(7)   \begin{equation*} \frac{\partial L}{\partial y} = 0 \end{equation*}

(7)式は、オイラー・ラグランジュ(Euler-Lagrange)方程式の特殊形である(一般形は、Ly'を含んでいるため、y'に関連する項が出る。変分ベイズにはy'は必要ない)。

次に条件2も同様に計算すると次のようになる。

(8)   \begin{eqnarray*} \delta I^{2} &=& \frac{1}{2} \epsilon^{2} \frac{d^{2} I}{d\epsilon^{2}} \biggm|_{\epsilon=0} \\ &=& \frac{1}{2} \epsilon^{2} \cdot \int \frac{d^{2} L}{d\epsilon^{2}} dx_{1}dx_{2} \biggm|_{\epsilon=0} \\ &=& \frac{1}{2} \int \frac{\partial^{2} L}{\partial y^{2}} (\epsilon \eta)^{2} dx_{1}dx_{2} \\ &<& 0 \end{eqnarray*}

(8)において、(\epsilon \eta)^{2} \ge 0なので第二変分が負となるためには、次の条件を満たせばよい。

(9)   \begin{equation*} \frac{\partial^{2} L}{\partial y^{2}} < 0 \end{equation*}

以上より、I[y]の変関数yに関する最大化問題は、(7),(9)の条件を満たす関数yを決定する問題に帰着する。なお、今回の導出では、yx_{1}, x_{2}の関数としたが、xがいくつあっても導出過程に代わりはない。

最後に、変関数yが次のような拘束条件を持つ場合についての解法を述べておく。

(10)   \begin{equation*} J = \int M(x_{1},x_{2},y) dx_{1}dx_{2} = K \end{equation*}

導出は同様なので省略するが、汎関数Iの代わりにI + \lambda (J - K)に対して、ラグランジュの未定乗数法で拘束条件を加えて変分問題を解くと、次の停留条件と極大値であるための条件を得る。

(11)   \begin{eqnarray*} G(x_{1},x_{2},y) &=& L(x_{1},x_{2},y) + \lambda M(x_{1},x_{2},y) \\ \frac{\partial G}{\partial y} &=& 0 \\ \frac{\partial^{2} G}{\partial y^{2}} &<& 0 \\ \end{eqnarray*}

\lambdaは拘束条件(10)から決定される。


以上で、変分法による汎関数最大化(極大化)問題の解き方が分かりました。次回は、この結果を利用して、変分ベイズの汎関数最大化問題の計算を行います。

変分ベイズの定式化から近似事後分布の導出まで(1)

PRML10.1の所で変分法を知らずに「なぜ?」と結構苦労したので、その備忘録です。

話の流れは以下になります。

  1. 変分ベイズの定式化(この記事)
  2. 準備:変分法
  3. 変分事後分布の導出

今日は1についてです。


1. 変分ベイズの定式化

注:全ての変数はスカラー、ベクトル、行列の何れかだが表記上は区別しない。また、離散変数の場合にも積和ではなく積分で一括表記する。

データDと潜在変数, パラメータ等のUnknownな全ての変数をZ=[Z_{1}, Z_{2}, ..., Z_{L}]とする。また、Zに関する任意の確率分布q(Z)を考える。この時、Dの対数周辺尤度について次の2つの関係が成り立つ。

(1)   \begin{equation*} \ln p(D) = \int q(Z) \ln p(D) dZ \end{equation*}

(2)   \begin{equation*} \ln p(D) = \ln p(D,Z) - \ln p(Z|D) \end{equation*}

(1)は1を掛けただけ、(2)はベイズの定理で対数をとっただけである。ここで、(2)を(1)に代入する。

(3)   \begin{eqnarray*} \lefteqn{ \ln p(D)  = \int q(Z) \{  \ln p(D,Z) - \ln p(Z|D)  \} dZ } \\ &=& \int q(Z) \{  \ln p(D, Z) - \ln p(Z|D) + \ln q(Z) - \ln q(Z)  \} dZ \\ &=& \int q(Z) \frac{ \ln p(D, Z) }{ \ln q(Z) } dZ + \int q(Z) \frac{ \ln q(Z) }{ \ln p(Z|D) } dZ \\ &=& F[q(Z)] + KL(q(Z) \| p(Z|D)) \\ \lefteqn{ F[q(Z)] =  \int q(Z) \frac{ \ln p(D, Z) }{ \ln q(Z) } dZ } \end{eqnarray*}

第二項KLZについての任意の分布q(Z)と真の事後分布p(Z|D)のカルバック・ライブラー情報量であり、2つの分布の近さを表す。KL \ge 0で二つの分布が一致する時にのみ0となる。

真の事後分布は未知なのでKLは操作できないが、(3)で定義した汎関数Fq(Z)に関して最大化することで、q(Z)p(Z|D)へとカルバック・ライブラー情報量の観点から接近し、結果として得られるq(Z)は事後分布の近似となる。よって、変分ベイズとは次の問題として定式化される。

(4)   \begin{equation*} p(Z|D) \approx \arg \max F[q(Z)] \end{equation*}

20130119a_VB

Fは変分下界と呼ばれ、対数周辺尤度ln p(D)の下限値を与えている。また、q(Z)のことを変分事後分布と呼ぶ。定式化は以上だが少し補足する。ここで、Fを分解してみる。

(5)   \begin{eqnarray*} F[q(Z)] &=&  \int q(Z) \frac{ \ln p(D, Z) }{ \ln q(Z) } dZ \\ &=&  \int q(Z) \ln p(D, Z) dZ + \int q(Z) \frac{ 1 }{ \ln q(Z) } dZ \\ &=&  \int q(Z) \ln p(D, Z) dZ - \int q(Z) \ln q(Z) dZ \\ \end{eqnarray*}

第1項の\ln p(D,Z)はUnknownな変数が与えられた時の対数尤度で、完全データの対数尤度と呼ぶ。完全データとは[D,Z]を指している。第1項の最大化は、「データとUnknownな変数の尤もらしい組合せを、期待値として評価して探せ」ということになる。一方で、第2項は変分事後分布q(Z)の平均エントロピー、すなわち分布のフラットさを表している(離散変数であれば、一様分布が平均エントロピー最大)。よって、第2項の最大化は、「変分事後分布が出来るだけフラットになるようにせよ」ということになる。

完全データの対数尤度の期待値は、変分事後分布がシャープなほど大きくなるが、エントロピーは変分事後分布がシャープなほど小さくなる。第1項と第2項は互いに対立する評価値で、変分ベイズはバランス点を探していることになる。

なぜ第1項だけを考えてはいけないか?それが変分ベイズ導入のモチベーションである。第一項を最大化するZを一点で推定するアルゴリズムをEMアルゴリズムと呼び、最尤法によるオーバーフィッティング問題が顕在化するためである。変分ベイズにおいては第2項が、オーバーフィッティングを回避するための正則化の役割を果たしている。


以上が変分ベイズによる事後分布近似方法の定式化になります。この解き方はまた今度。