はじめに
以前、ベイズ推論による機械学習入門という書籍を読み、以下のような記事を公開しました。
この書籍には、多次元ガウス混合モデルによるクラスタリング手法が詳しく述べられていますが、多次元ガウス混合モデルによる教師あり分類の手法は少し触れる程度にとどまっています(クラスタリングで述べられている導出過程を流用すれば、比較的簡単に導出できるためだと思います)。
ということで、今回は多次元ガウス混合モデルによる教師あり分類の導出、実装をやってみました。
ガウス混合モデル
モデル自体は教師ありだろうが、教師なしだろうが同じモデルを考えることができます。まず、全て独立でサンプリングされる次元のデータがと個あるとし、そのそれぞれに対応するラベルデータをとします。ここで、は of 表現と呼ばれ、クラス数がの場合、所属するクラスに対応する要素のみで他は全ての次元ベクトルです。このをサンプルするための分布として、パラメータを持つカテゴリ分布を考えます。
また、個のガウス分布が混合したモデルを考えるので、それぞれの分布の平均を、精度行列(共分散行列の逆行列)をとします。なお、各パラメータを個別に扱わない場面では表記を簡単にするため、これらのガウス混合分布のパラメータを全てまとめてと表記します。
更に、カテゴリ分布のパラメータや各ガウス分布の、といったパラメータもベイズ推論の枠組みでは、確率分布に従うと考えます。ここでは、それぞれ共役事前分布に従うことにします。具体的には、
という分布に従うとします。ここで、はディリクレ分布、はガウス・ウィシャート分布です(具体的な分布の式は書籍を参照ください)。ここで事前分布として共役事前分布を採用しておくと事後分布や予測分布が解析的に求めることができます。やは超パラメータで今回は事前に与えるものとします。また、については、それぞれのガウス分布で異なる値を与えることもできますが、簡単のため今回は全て共通の値とします。
さて、ここまで出てきた確率変数の同時分布を考えると、以下のようになります。
ここで、はガウス分布、はカテゴリ分布を表し、とは式(1)、式(2)で表される分布です。
このモデルのグラフィカルモデルを下図に示します。この図でとの四角で囲まれているのが学習データとして与えられるデータ、とがクラスを予測する対象のデータ(テストデータ)です。
このように図示しておくと有効分離の考え方を用いて、変数間の条件付き独立性がわかりやすくなり、後々式展開する際に便利です。有効分離については、以下の記事がわかりやすいと思います。
推論
上述式(3)のモデルから、パラメータの事後分布とそれを元に、新たなデータが与えられた時のの確率分布である予測分布を導出します。このモデルの場合、上述のように共役事前分布を事前分布として採用しているため、事後分布および予測分布は解析的に求めることができます。ここでは、逐次的にデータが与えられるとして、オンラインで予測分布を更新していけるような式を導出します。
事後分布の導出
今回扱う問題ではとがデータとして与えられるので、これらが与えられている下でのパラメータの事後分布を求めます。つまり、
を求めることになります。なお、クラスタリングの場合、がデータとして与えられているわけではないので、求めるべき事後分布としてはとなり、こちらは解析的に求めることが難しくなります(なので、ギブスサンプリングや変分推論などの近似手法を使う)。
ここで式(3)を見ると、とは完全に分離していることがわかります。つまり、とが与えられた下では、とは独立であり、
を求めればよいことになります。
それでは、式(5)を展開してみます。式(5)の右辺について、対数をとってみると、
となります。これを見ると各に対して分離しているのがわかるので、各に対しそれぞれ独立に事後分布を求めればよいことになります。また、式(7)の中括弧の中はガウス・ウィシャート分布を事前分布として用いたガウス分布パラメータの事後分布推定の計算となります。この計算方法は件の書籍に詳しく載っているので、ここでは割愛しますが、最終的な形は以下のようになります。
ただし、
です。同様に式(6)はディリクレ分布を事前分布としてカテゴリ分布の事後分布そのものであるため、
というディリクレ分布になります。ただし、
です。
ここで、式(9)および式(11)の各超パラメータをデータが与えられるたびに逐次更新できる形に変更すると、
となり、これによりの時点で、計算した超パラメータを用いて、番目のデータが来たときの超パラメータに更新することが可能です(初期値は事前分布の超パラメータ)。
予測分布の導出
最終的に欲しいのは、新たなデータが与えられた時のクラス分類結果です。そのため、が与えられた下でのの予測分布を求めます。予測分布は、
となります。ここで、に関係のない項は無視していいので、とは考慮から外します。更に、について、グラフィカルモデルから有効分離性を読み解くとが与えられた下では、とは独立であり(に繋がるが観測され、かつ、head-to-head型を形成しているにおいてが観測されていないため)、となるため、こちらも考慮しなくてよくなります。まとめると、式(13)は
と、2つの項のみ考えればよくなります。1つ目の項は、
であり、ここでもグラフィカルモデルから有効分離性を考えると、とすることができます。従って、
となります。
は、式(5)、式(8)、式(9)に記している事後分布そのものです。今、各に対してであるそれぞれの予測分布を求めたいので、式(16)より
を求めればよいことになります。これは、ガウス・ウィッシャート分布を事前分布として用いた場合のガウス分布の予測分布計算になります。こちらも予測分布導出は書籍に載っているので、導出過程は省略し最終形を示します。
ただし、は多次元のスチューデントのt分布で、はの次元です。
次に、式(14)のですが、こちらも同様に、
となります。ただし、第2式から第3式への展開はグラフィカルモデルよりを観測した下でのとの条件付き独立性を利用しました。そしてこの式(19)は、ディリクレ分布を事前分布として使ったカテゴリ分布の予測分布です。ということで、この予測分布についても導出過程は割愛し最終形を示します。
ただし、の各要素は
です。
ここまでの話をまとめると、番目のデータが与えら手た時の予測分布の更新式は、式(14)、式(17)、式(20)、式(21)より
となります(各超パラメータは式(12)、式(21))。そして、各クラスの所属確率の平均を得たいときは、に対し、式(22)の右辺にを代入し、得られる値を合計が1になるように正規化すればよいです。
実装
ということで実装です。
import numpy as np import scipy as sci import math import random import matplotlib.pyplot as plt import matplotlib.animation as anm class multivariate_student_t(): def __init__(self, mu, lam, nu): self.D = mu.shape[0] self.mu = mu self.lam = lam self.nu = nu def pdf(self, x): temp1 = np.exp( math.lgamma((self.nu+self.D)/2) - math.lgamma(self.nu/2) ) temp2 = np.sqrt(np.linalg.det(self.lam)) / (np.pi*self.nu)**(self.D/2) temp3 = 1 + np.dot(np.dot((x-self.mu), self.lam), (x-self.mu))/self.nu temp4 = -(self.nu+self.D)/2 return temp1*temp2*((temp3)**temp4) class mixture_gaussian_classifier: def __init__(self, alpha0, m0, invW0, beta0, nu0): self.K = alpha0.shape[0] #予測分布(多次元t分布)の生成 self.pred_distr = (lambda m, W, beta, nu: multivariate_student_t (m, (1 - m.shape[0] + nu) * beta * W/(1 + beta), (1 - m.shape[0] + nu))) #超パラメータの初期化 self.alpha = alpha0 self.eta = alpha0/np.sum(alpha0) self.m = [None] * self.K self.invW = [None] * self.K self.beta = np.zeros(self.K) self.nu = np.zeros(self.K) for k in range(self.K): self.m[k] = m0 self.invW[k] = invW0 self.beta[k] = beta0 self.nu[k] = nu0 def train(self, x, s): #各超パラメータの更新 k = np.where(s == 1)[0][0] self.alpha[k] += 1.0 self.eta = self.alpha/np.sum(self.alpha) pre_beta = self.beta[k] self.beta[k] += 1.0 self.nu[k] += 1.0 pre_m = self.m[k] self.m[k] = (x + pre_m * pre_beta)/self.beta[k] self.invW[k] = (self.invW[k] + np.dot(x.reshape(-1,1), x.reshape(-1,1).T) - self.beta[k] * np.dot(self.m[k].reshape(-1,1), self.m[k].reshape(-1,1).T) + pre_beta * np.dot(pre_m.reshape(-1,1), pre_m.reshape(-1,1).T)) def test(self, x): s = np.zeros(self.K) #各クラスの予測平均を算出 for k in range(self.K): s[k] = (self.eta[k] * self.pred_distr(self.m[k], np.linalg.inv(self.invW[k]), self.beta[k], self.nu[k]).pdf(x)) s /= np.sum(s) return s
クラスmixture_gaussian_classifierは、ここまで述べてきた学習(train)と予測平均の算出(test)を行うクラスです。また、multivariate_student_t多次元のスチューデントのt分布のクラスで、この実装については、以下の記事の内容を参考にさせてもらいました。
実験・結果
ということで、以下のようなコードを組んで実験してみました。クラス数はで、学習データ数はそれぞれと少し偏りを持たせています。
また、せっかく事後分布を逐次学習できるようにしたので、学習データが与えられるたびに予測平均がどのように変化するかをアニメーションにしてみました。
global used_idx def update(i, idx, mgc, X, S): plt.cla test_range = range(0,60) X_new = np.zeros((2, len(test_range)**2)) S_new = np.zeros((3, len(test_range)**2)) i = idx[i] #データの学習 mgc.train(X[:, i],S[:, i]) used_idx.append(i) c = S[:, used_idx] plt.scatter(X[0, used_idx], X[1, used_idx], edgecolors="black", marker="o" , color=S[:, used_idx].T) count = 0 S_res = np.zeros((len(test_range), len(test_range), 3)) for x1 in test_range: for x2 in test_range: #各座標位置に対する予測平均を算出 X_new[:, count] = np.array([float(x1), float(x2)]) S_new[:, count] = mgc.test(X_new[:, count]) S_res[int(x2 - test_range[0]), int(x1 - test_range[0]), :] = S_new[:, count] count += 1 plt.imshow(S_res) def main(): N1 = 50 N2 = 60 N3 = 70 X = np.zeros((2, N1+N2+N3)) S = np.zeros((3, N1+N2+N3), 'int8') #クラス1の学習データ mu1 = np.array([45, 25]) A = np.array([[3,0],[0,8]]) sigma1 = np.dot(A.T, A) X[:, :N1] = np.random.multivariate_normal(mu1, sigma1, N1).T S[0,:N1] =1 #クラス2の学習データ mu2 = np.array([35, 40]) A = np.array([[5, 3],[3,7]]) sigma2 = np.dot(A.T, A) X[:, N1:N1 + N2] = np.random.multivariate_normal(mu2, sigma2, N2).T S[1,N1:N1 + N2] = 1 #クラス3の学習データ mu3 = np.array([30, 25]) A = np.array([[7,-5],[-5,7]]) sigma3 = np.dot(A.T, A) X[:, N1 + N2:]= np.random.multivariate_normal(mu3, sigma3, N3).T S[2,N1 + N2:] = 1 #事前分布の超パラメータ alpha0 = np.random.randn(3) m0 = (mu1 + mu2 + mu3)/2 invW0 = np.linalg.inv((sigma1 + sigma2+ sigma3)/3) beta0 = 1.0 nu0 = m0.shape[0] #推論オブジェクトの生成 mgc = mixture_gaussian_classifier(alpha0, m0,invW0, beta0, nu0) fig = plt.figure() idx = [i for i in range(X.shape[1])] random.shuffle(idx) global used_idx used_idx = [] #結果のアニメーション表示(300msごとに関数updateを呼び出す) ani = anm.FuncAnimation(fig, update, fargs = (idx, mgc, X, S), interval = 300, frames = N1 + N2+ N3) ani.save("result.gif", writer="imagemagick") if __name__ == '__main__': main()
そして、こちらが結果です。
この図で、各色がそれぞれのクラスを表し、追加されていく各点が学習データを表しています。また、背景色が予測分布が算出した各座標点の各クラスへの所属確率を表しています。
特に左上に辺りの領域に着目すると、最初の方は青色のクラスと予測されていましたが、青と緑のクラスのデータを観測していくにつれ、緑色のクラスだと予測されるようになっていくのがわかります。これは青色のクラスが左下から右上にかけての共分散が大きいのに対し、緑色のクラスの共分散が左上から右下への共分散が大きく、その特徴をサンプルが増えるにつれ、正しく推論できるようになっていくためだと考えられます。