甲斐性のない男が機械学習とか最適化とかを綴るブログ

うどんくらいしか食べる気がしない

MCMCでガウス過程分類(ロジスティック回帰編)

2020/2/2追記

ソースコードに大きな間違いを見つけたので修正した。2か所間違いがあってそれらが奇跡的に作用することでそれっぽい結果に見えていた。詳細は実装の節を参照。

はじめに

以下の記事の中で、ガウス過程と機械学習サポートページにあるサンプルコードをもとにコーシー分布を観測モデルに基づくガウス過程回帰を実装(というより微修正)した。

yamagensakam.hatenablog.com

今回はこれを少しだけ変えて、ガウス過程ロジスティック回帰をMCMCで行い2クラス分類問題を解いてみる。

やったこと

といっても、本質的には観測モデルp(y|f({\bf x}))

\displaystyle p(y|f({\bf x})) = {\rm Bern}(y| \sigma(f({\bf x}) ) = \left\{ \sigma(f({\bf x}))\right\}^y \left\{1 - \sigma(f({\bf x}))\right\}^{(1 - y)}

というベルヌーイ分布にするだけ。あとはテスト入力に対するyの予測値は、テスト入力{\bar {\bf x}}に対する関数出力値の事後分布からのサンプリング値f_1({\bar {\bf x}}), \cdots,f_S({\bar {\bf x}})を用いて、

\displaystyle{\bar y} = \frac{1}{S} \sum_{i = 1}^{S}  \sigma(f_i({\bar {\bf x}}))

として近似する。

実装

以下、ソースコード

2020/2/20追記

冒頭でも書いたが大ポカを犯していた。間違いは

  • ベルヌーイ分布の対数尤度にマイナスをつけていた。
  • 学習データをプロットするときに、クラスの色を逆にしていた。

という2点。対数尤度にマイナスをつけてしまっていたので予測値としては逆の結果が得られるが、プロットする際クラスラベルを逆にするというミスも犯していたから結果的にそれっぽい出力が得られ気が付かなかった。ということで、修正したソースコードは以下の通り。

import sys
import numpy as np
from pylab import *
from numpy import exp,log
from elliptical import elliptical
from numpy.linalg import cholesky as chol

def gpr_bell(f, param):
    X,y,Kinv = param
    N = X.shape[1]
    return gpr_bell_lik (y[0:N], f[0:N], Kinv)

def gpr_bell_lik(y, f, Kinv):
    log_sigf = log(1.0/(1.0+exp(-f)))

    # 2020/2/2 ベルヌーイ分布の対数尤度にマイナスをつけていたがない方が正しいので修正
    # return -(np.sum(y @ log_sigf) + np.sum((1 - y) @ (1 - log_sigf)))\
    #         - np.dot (f, np.dot(Kinv, f)) / 2            
    return np.sum((y * log_sigf) + ((1 - y) * (1 - log_sigf)))\
            - (f @ Kinv @ f / 2)

def GetDistanceMatrix(X):
    H = np.tile(np.diag((X.T @ X)) , (np.size(X , 1) , 1))
    G = X.T @ X
    DistMat = H - 2 * G + H.T
    return DistMat

def kgauss(tau, sigma):
    return lambda X: tau * exp (-GetDistanceMatrix(X) / sigma)

def kernel_matrix(X, kernel):
    N = X.shape[1]
    eta = 1e-6
    return kernel(X) + eta * np.eye(N)


def gpr_mcmc(X, y, cov_func, iters, remove_num):
    
    N = X.shape[1]
    K_NN = kernel_matrix(X, cov_func)
    K_NNinv = inv(K_NN)

    #こうすると多次元正規分布からサンプリングすることになる
    P = chol (K_NN)
    f = P @ randn(N)
    
    #MCMCサンプルのうち最初のremove_num個は使わない
    S = iters - remove_num
    f_posterior = np.zeros((f.shape[0], S))

    for iter in range(iters):
        #楕円スライスサンプリング
        f,lik = elliptical (f, P, gpr_bell, (X, y, K_NNinv))
        print ('\r[iter %2d]' % (iter + 1),)
        if iter >= remove_num:
            f_posterior[:, iter - remove_num] = f
    return K_NNinv, f_posterior

def predict(X_train, X_test, cov_func, K_NNinv, f_posterior, sampling_num):
    N = X_train.shape[1]
    M = X_test.shape[1]
    XX = np.concatenate([X_train, X_test], axis = 1)

    K = kernel_matrix(XX, cov_func)
    K_NM = K[:N, N:N + M]
    K_MM = K[N:N + M, N:N + M]
    y_test = np.zeros(M)
    S = f_posterior.shape[1]

    for t in range(sampling_num):
        s = np.random.randint(S)
        f_n = f_posterior[:, s]
        mu = K_NM.T@ K_NNinv @ f_n
        sig = K_MM - K_NM.T @ K_NNinv @ K_NM
        P = chol(sig)
        f_m = P @ np.random.randn(M) + mu
        y_test += (1.0/(1.0+np.exp(-f_m)))

    y_test /= sampling_num
    return y_test

def make_data(xmin, xmax, n):
    np.random.seed(2222)
    cnt = 0
    t = 0
    N = n * 9
    X_train = np.zeros((2, N))
    y_train = np.zeros(N)
    center1 =  np.linspace(-(-xmin[0] + xmax[0])/3.5, (-xmin[0] + xmax[0])/3.5, 3)
    center2 =  np.linspace(-(-xmin[1] + xmax[1])/3.5, (-xmin[1] + xmax[1])/3.5, 3)
    for i in range(3):
        for j in range(3):
            y_train[t:t + n] = cnt % 2
            X_train[:, t:t + n] = 1.5 * np.random.randn(2, n) \
                                    + np.array([center1[i], center2[j]])[:, None]
            t += n
            cnt += 1


    #テスト入力点生成
    T1 = np.linspace(xmin[0], xmax[0], 60)
    T2 = np.linspace(xmin[1], xmax[1], 60)
    X_test = np.zeros((2, T1.shape[0] * T2.shape[0]))
    i = 0
    for t1 in T1:
        for t2 in T2:
            X_test[:, i] = np.array([t1, t2])
            i += 1
    return X_train, y_train, X_test
    
def main ():    
    xmin = np.array([-9, -9])
    xmax = np.array([9, 9])
    
    #学習データとテストデータの生成
    n = 20
    X_train, y_train, X_test = make_data(xmin, xmax, n)

    #共分散関数の設定
    cov_func = kgauss(3, 2.5)
    
    #MCMCで事後分布からのサンプリング
    iters = 1000
    remove_num = 200
    K_NNinv, f_posterior = gpr_mcmc(X_train, y_train, cov_func, iters, remove_num)

    #予測値の計算
    sampling_num = 50
    y_test = predict(X_train, X_test, cov_func, K_NNinv, f_posterior, sampling_num)

    #結果の表示
    size = 100
    sc = plt.scatter(X_test[0, :], X_test[1, :], size, y_test, marker = ',', cmap='bwr')
    plt.colorbar(sc)

    # 2020/2/2 プロットの色が逆だったので修正
    # plt.scatter(X_train[0, y_train == 1], X_train[1, y_train == 1] ,
    #             edgecolors="black", marker="o" , color='b')
    # plt.scatter(X_train[0, y_train == 0], X_train[1, y_train == 0] ,
    #             edgecolors="black", marker="o" , color='r')
    plt.scatter(X_train[0, y_train == 1], X_train[1, y_train == 1] ,
                edgecolors="black", marker="o" , color='r')
    plt.scatter(X_train[0, y_train == 0], X_train[1, y_train == 0] ,
                edgecolors="black", marker="o" , color='b')

    plt.axis([xmin[0], xmax[0], xmin[1], xmax[1]])

    plt.show()

if __name__ == "__main__":
    main ()

上で書いたこと以上に前回のコードから変更しているが、本質的な変更は尤度計算がgpr_bell_likになっておりベルヌーイ分布を使った計算になっていること。他は見れば大体わかると思うが、1つだけ補足するとカーネル行列を作るところで以前はfor文を回していたが遅いのでfor文使わないように修正している。詳細はこの記事の『各データ間の距離の2乗を計算』ってとこ参照。

yamagensakam.hatenablog.com

結果

f:id:YamagenSakam:20200202021523p:plain

丸マーカーでプロットしているのが学習データ、これをもとに各座標位置の所属クラスを予測している。ちゃんと非線形なクラス分類ができている(2020/2/2 修正に伴い当然結果も変わる)。

最後に

ロジスティック回帰編と書いているってことは、プロビット回帰編もいつかやる・・・はず。

ガウス過程回帰で学習データにない任意の入力点に対する関数出力値をMCMCサンプリングする方法

はじめに

表題についてちょっと調べた&やってみたのでメモ。 やりたいことは表題の通りなんだけど、少し数式で表現して整理。

  • 今、学習データとして \left\{{\bf X}, {\bf y} \right\} = \left\{({\bf x}_1, y_1), \cdots, ({\bf x}_N, y_N)\right\}N個得られているとする。
  • 入力{\bf x}と関数fが与えられたとき、出力yp(y | f({\bf x}))に従うとする(ハイパーパラメータを持つ場合もあるが今回は固定値として扱うので表記から省略)。
  • そして、この関数fガウス過程p(f)={\mathcal gp}(0, k({\bf x}, {\bf x}'))に従う確率変数とする( k({\bf x}, {\bf x}')は共分散関数)。
    • p(y | f({\bf x}))ガウス分布でないなら、関数の事後分布p(f|{\bf X}, {\bf y}) \propto \prod_{i = 1}^{N}p( y_i | f({\bf x}_i))p(f)は解析的に得られないのでそういう場合は近似するかMCMC等でサンプリングを行う。
  • やりたいのは、M個の任意の未知入力\left\{{\bar {\bf x}}_1, \cdots, {\bar{\bf x}}_M \right\}が与えられて、事後分布からサンプリングされた関数fの出力\left\{f({\bar {\bf x}}_1), \cdots, f({\bar{\bf x}}_M) \right\}を得たい。
    • 予測分布からのサンプリングも得たいという話もあるけど今回は対象外。
      • とりあえず、事後分布からの出力値f({\bar {\bf x}})が得られれば、これからp(y |f({\bar {\bf x}}))が計算できるので、それっぽいyは得られそう。

ようするに、学習入力点N個に対する関数fの出力を並べたベクトルを{\bf f}_N = \left[f({\bf x})_1 \cdots, f({\bf x}_N)\ \right]^Tとし、未知入力点M個に対する関数fの出力を並べたベクトルを{\bar {\bf f}}_M = \left[f({\bar {\bf x}})_1 \cdots, f({\bar {\bf x}}_M)\ \right]^Tとした場合、

\displaystyle p(f|{\bf X}, {\bf y})  =  p({\bf f}_N | {\bf y}) \propto p({\bf y}| {\bf f}_N) p({\bf f}_N)

からのMCMCサンプリングはできるのは分かるけど、 p({\bar {\bf f}}_M |  {\bf y}) からのサンプリングはどうやったらいいんだ?って話。

方法1 ~MCMCサンプリングの時、同時分布 p({\bf f}_N, {\bar {\bf f}}_M )からサンプリング~

非常にお世話になっている書籍

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

サポートページに、コーシー分布を観測モデルとしたガウス過程回帰のpythonコードgpr-cauchy.py, elliptical.pyがある。これは楕円スライスサンプリングというMCMCサンプリングの一手法を用いて関数出力値の推定を行っているが、これを見るとサンプルを得るところで


p({\bf f}_N, {\bar {\bf f}}_M ) = {\mathcal N} \left( {\bf f}_N, {\bar {\bf f}}_M| 0, \begin{pmatrix}
{\bf K}_{NN} & {\bf K}_{NM}  \\
{\bf K}_{NM}^T & {\bf K}_{MM} \\
\end{pmatrix}
\right)

(ただし、{\bf K}_{NN}K_{i,j} = , k({\bf x}_i, {\bf x}_j) \; (i,j \in \left\{1, \cdots, N \right\})という要素を持つN \times Nの行列、{\bf K}_{NM}K_{i,j} = , k({\bf x}_i, {\bar {\bf x}}_j) \; (i \in \left\{1, \cdots, N \right\}, j \in  \left\{1, \cdots, M \right\})という要素を持つN \times Mの行列、{\bf K}_{MM}K_{i,j} = , k({\bar{\bf x}}_i, {\bar{\bf x}}_j) \; (i,j \in \left\{1, \cdots, M \right\})という要素を持つM \times Mの行列)

からのサンプルを使って更新している。そして、そのサンプルを採択するか否かの尤度計算はp({\bf y}| {\bf f}_N) p({\bf f}_N)で計算、つまり学習入力点だけで行っている*1

なるほど、確かにこうすればp({\bf f}_N, {\bar {\bf f}}_M |  {\bf y}) からサンプリングしたことになり、そのサンプルから{\bf f}_Nを捨てれば周辺化することになるからp({\bar {\bf f}}_M |  {\bf y}) からサンプリングしたことになる。

方法2 ~ p({\bf f}_N | {\bf y}) からのサンプルを用いて p({\bar {\bf f}}_M |  {\bf y}) からサンプリング~

方法1だと、学習時に評価したい未知入力\left\{{\bar {\bf x}}_1, \cdots, {\bar{\bf x}}_M \right\}がわかっている必要があるが、場合によってはそれが困難だったりする。そういう場合はどうするか?

p({\bf f}_N | {\bf y})からのサンプルがS{\bf f}_N^{(1)}, \cdots, {\bf f}_N^{(S)}得られているとすると、

\displaystyle p({\bar {\bf f}}_M |  {\bf y}) = \int p({\bar {\bf f}}_M | {\bf f}_N) p({\bf f}_N | {\bf y})d  {\bf f} _N  \sim \frac{1}{S} \sum_{s = 1}^{S}p({\bar {\bf f}}_M | {\bf f}_N^{(s)})

と近似できる。ここでガウス過程の性質より、

p({\bar {\bf f}}_M | {\bf f}_N^{(s)}) = {\mathcal N} ({\bar {\bf f}}_M | {\bf K}_{NM}^T  {\bf K}_{NN}^{-1}  {\bf f}_N^{(s)}  , {\bf K}_{MM} -  {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf K}_{NM} )

なので、p({\bar {\bf f}}_M |  {\bf y})は全て同じ混合比の混合正規分布で近似できる。従ってs \in \left\{1 \cdots S \right\}の乱数を発生させた後、それに紐づく多次元正規分布 {\mathcal N} ({\bar {\bf f}}_M | {\bf K}_{NM}^T  {\bf K}_{NN}^{-1}  {\bf f}_N^{(s)}  , {\bf K}_{MM} -  {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf K}_{NM} )からサンプリングすればよい。

実装

方法2によるコーシー分布ガウス過程回帰を上記サポートページにあるgpr-cauchy.pyを改変する形でやってみた。

import sys
import putil
import numpy as np
from pylab import *
from numpy import exp,log
from elliptical import elliptical
from numpy.linalg import cholesky as chol

def gpr_cauchy (f, param):
    x,y,gamma,Kinv = param
    M = len(x)
    return gpr_cauchy_lik (y[0:M], f[0:M], gamma, Kinv)

def gpr_cauchy_lik (y, f, gamma, Kinv):
    return - np.sum (log (gamma + (y - f)**2 / gamma)) \
           - np.dot (f, np.dot(Kinv, f)) / 2

def kgauss (tau,sigma):
    return lambda x,y: exp(tau) * exp (-(x - y)**2 / exp(sigma))

def kernel_matrix (xx, kernel):
    N = len(xx)
    eta = 1e-6
    return np.array (
        [kernel (xi, xj) for xi in xx for xj in xx]
    ).reshape(N,N) + eta * np.eye(N)

def gpr_mcmc(x, y, iters, xmin, xmax, gamma):
    
    N = len(x)
    K_NN = kernel_matrix(x, kgauss(1,1))
    K_NNinv = inv(K_NN)

    #こうすると多次元正規分布からサンプリングすることになる
    P = chol (K_NN)
    f = P @ randn(N)
    
    #MCMCサンプルのうち最初の200個は使わない
    remove_num = 200
    S = iters - remove_num
    f_posterior = np.zeros((f.shape[0], S))
    
    for iter in range(iters):
        #楕円スライスサンプリング
        f,lik = elliptical (f, P, gpr_cauchy, (x,y,gamma,K_NNinv))
        print ('\r[iter %2d]' % (iter + 1),)
        if iter >= remove_num:
            f_posterior[:, iter - remove_num] = f
    
    #未知入力点
    x_new = np.linspace(xmin, xmax, 100)
    M = len(x_new)
    xx = np.hstack((x, x_new))

    K = kernel_matrix(xx, kgauss(1, 1))
    K_NM = K[:N, N:N + M]
    K_MM = K[N:N + M, N:N + M]

    g = np.zeros (len(x_new))
    #50個サンプリング
    sampling_num = 50
    for t in range(sampling_num):
        s = np.random.randint(S)
        f_n = f_posterior[:, s]
        mu = K_NM.T@ K_NNinv @ f_n
        sig = K_MM - K_NM.T @ K_NNinv @ K_NM
        P = chol(sig)
        f_m = P @ np.random.randn(M) + mu
        g = g + f_m
        plot (x_new, f_m)
    plot (x, y, 'bx', markersize=14)
    plot (x_new, g/sampling_num, 'k', linewidth=3)

def usage ():
    print ('usage: gpr-cauchy.py data.xyf iters [output]')
    sys.exit (0)

def main ():
    xmin = -5
    xmax =  5
    ymin = -7.5
    ymax = 12.5
    gamma = 0.2
    
    if len(sys.argv) < 3:
        usage ()
    else:
        [x, y, f] = np.loadtxt (sys.argv[1]).T
        iters = int (sys.argv[2])

    gpr_mcmc (x, y, iters, xmin, xmax, gamma)
    axis ([xmin, xmax, ymin, ymax])
    show ()

if __name__ == "__main__":
    main ()

楕円スライスサンプリングを行う関数ellipticalについては、 サポートページのelliptical.pyから変えてないので省略。また多次元正規分布からはコレスキー分解を用いてサンプリングしている(この辺が参考になる)。これを、

python .\gpr-cauchy_new.py  .\gpr-cauchy.dat 1000

として、実行すると以下のような結果が得られた。

f:id:YamagenSakam:20200123235146p:plain

うん、いい感じ。

今後やりたいこと

  • ハイパーパラメータも事前分布与えてMCMCサンプリング
  • 予測分布からのサンプリング
  • ガウス過程分類
  • 変分近似

*1:事前分布の方はp({\bf f}_N, {\bar {\bf f}}_M)の同時分布を評価するべきかと思ったがどうなんだろ?実際両方やってみたけど、結果はあまり変わらなかった

EMアルゴリズムと変分ベイズEMアルゴリズムのメモ

表題の通り。この辺りややこしくていつもわからなくなるので、身につけるためにも自分なりの理解をまとめたメモを書く。

対数尤度と下限とKLダイバージェンス

まず、既知の変数をX、未知の変数(ここでは潜在変数か未知パラメータかは区別していない)をZとすると以下が成り立つ。

\displaystyle \ln p(X) = KL\left[q(Z)\parallel p(Z| X) \right] + {\mathcal L}\left[q(Z)\right]  \tag{1}

ただし、KL\left[q(Z) \parallel p(Z | X) \right] はKLダイバージェンスで、

\displaystyle KL\left[q(Z) \parallel  p(Z | X) \right] =\int q(Z) \ln \frac{q(Z)}{p(Z | X)} dZ = - \int q(Z) \ln \frac{p(Z | X)}{q(Z)} dZ \tag{2}

{\mathcal L}\left[q(Z)\right] は、対数尤度の下限で

\displaystyle {\mathcal L}\left[q(Z)\right] = \int q(Z) \ln \frac{p(X, Z)}{q(Z)} dZ \tag{3}

であり、常に

\displaystyle \ln p(X) \geq {\mathcal L}\left[q(Z)\right] \tag{4}

が成り立つ。これは関数\lnが上に凸でありイェンセン不等式より\int p(x)\ln f(x)dx \leq \ln \int p(x)f(x)dxであるため、\int q(Z) \ln \frac{p(X, Z)}{q(Z)} dZ \leq \ln \int q(Z)\frac{p(X, Z)}{q(Z)} dZ =\ln \int p(X, Z) dZ =  \ln p(X) として得られる。

また、式(1)は、式(2)と式(3)を実際に足してみると得られる。

EMアルゴリズム

考え方

  • 得られているものは既知の変数のX、求めたいのは未知パラメータ\Theta
  • 目的関数は対数尤度\ln p(X|\Theta)であり、これを最大化する\Thetaを求めたい。
    • ところが、これを直接的に最適化するのは一般に困難。
  • ここで、潜在変数Zを導入し、{\mathcal L}\left[q(Z), \Theta \right] = \int q(Z) \ln \frac{p(X, Z|\Theta)}{q(Z)} dZという量を考える。
    • 式(4)より、\ln p(X|\Theta) \geq {\mathcal L}\left[q(Z), \Theta \right]が成り立つ。
  • そのため、{\mathcal L}\left[q(Z), \Theta \right]を最大化する関数qとパラメータ\Thetaを求めることを考える。
    • ただし、この両者を一気に最適化する解析解は得られないので、座標降下法の考えで一方を固定し他方を最適化するという操作を繰り返して求める。

EMアルゴリズムの更新式

ということで、更新アルゴリズムは以下の通り。

  1. まずは、パラメータを既知の値 {\hat \Theta}と固定した場合に、{\mathcal L} \left[q(Z), {\hat \Theta} \right]を最大にする関数qを求める。
    • 式(1)より、 \displaystyle {\mathcal L}\left[q(Z), {\hat \Theta} \right]  = \ln p(X| {\hat \Theta}) - KL\left[q(Z) | p(Z| X, {\hat \Theta}) \right]
    • KL\left[q(Z) \parallel p(Z| X, {\hat \Theta})\right] \geq 0 なので、結局q(Z) = p(Z| X, {\hat \Theta})とするのが最適。
    • 色々分け方があるようだが、これをEステップと呼ぶことが多いらしい。
  2. 次に、関数qを上で求めたq(Z) = p(Z| X, {\hat \Theta})に固定し、{\mathcal L} \left[p(Z| X, {\hat \Theta}),  \Theta \right]を最大化する\Thetaを求める。
    • {\mathcal L} \left[p(Z| X, {\hat \Theta}),  \Theta \right] =  \int p(Z| X, {\hat \Theta})\ln \frac{p(X, Z|\Theta)}{p(Z| X, {\hat \Theta})} dZ =  \int p(Z| X, {\hat \Theta})\ln p(X, Z|\Theta) dZ + {\rm const}なので、 \int p(Z| X, {\hat \Theta})\ln p(X, Z|\Theta) dZ = \langle \ln p(X, Z|\Theta)\rangle _{p(Z| X, {\hat \Theta})}を最大化する\Thetaを求めればよい。
    • 期待値の最適化なので\ln p(X, Z|\Theta)さえ最適化できるなら適用可能。
    • こちらをMステップと呼ぶ。

変分ベイズEMアルゴリズム

変分近似の考え方

  • 得られているものは既知の変数のX、未知変数は潜在変数もパラメータもひっくるめてZとする。
  • モデルp(X, Z)が与えられていて、未知変数の事後分布p(Z | X)を求めたい。
    • これは、p(Z | X) = \frac{p(X, Z)}{p(X)}であり、厳密に求めるのは困難(p(X)が難しい)。
  • そこで、p(Z | X)を別の分布q(Z)で近似することを考える。
    • 当然q(Z)p(Z | X)に近い方がいいので、式(2)のKLダイバージェンスを最小化する関数qを求めることになる。
  • その際、q(Z)に何の制約も課さないと、q(Z) = p(Z | X) という無意味な解が得られるので、q(Z)に何らの制約を課して表現能力を限定する。
    • 典型的なのは、q(Z) = \prod_{i=1}^{D} q(z_i | X)のように未知変数間の独立性を制約にする近似(平均場近似)。
  • 式(1)より、KLダイバージェンスを最小化は{\mathcal L}\left[q(Z)\right] を最大化することと等価。
    • \ln p(X)は定数なので。
  • ということで、定めた制約のもと{\mathcal L}\left[q(Z)\right] を最大化するqを求めればよい。

ここまでは変分近似の考え方。おそらく、変分ベイズEMアルゴリズムと言うと、以下のようなもう少し限定されたアルゴリズムのことを言うと思う。

変分EMアルゴリズムの考え方

  • 今、未知変数がパラメータ\Thetaと潜在変数Zだとする。
  • これらの事後分布p(\Theta, Z | X)を求めたいが厳密に求めるのは困難なので、q(\Theta, Z)という分布で近似することを考える。
    • qの制約としては、q(\Theta, Z) = q(\Theta)q(Z)という独立性を制約とする。
    • モデルに応じてどういう独立性を仮定するのが適切なのかが変わってくる(\Thetaをいくつかの要素が独立だと仮定したり)。
  • 式(3)より、\displaystyle {\mathcal L}\left[q(\Theta)q(Z)\right] = \int \int q(\Theta)q(Z) \ln \frac{p(X|\Theta, Z)p(Z | \Theta)p(\Theta)}{q(\Theta)q(Z)} dZ d\Theta
    • この式はベイズの定理よりp(X, \Theta, Z)=p(X|\Theta, Z)p(Z | \Theta)p(\Theta)と展開している。
      • これはベイズの定理から導かれる展開で、モデル次第ではもっと色々展開できる可能性がある(p(Z|\Theta)=p(Z)とできたりとか)。
  • 後はこれを座標降下法の要領で、q(\Theta)を固定してq(Z)に関する最適化とq(Z)を固定してq(\Theta)に関する最適化を交互に繰り返せばよい。

変分EMアルゴリズムの更新式の導出

ここからはもう少し、この式を展開する。


\begin{align}
{\mathcal L} \left[q(\Theta)q(Z) \right] &= \int \int q(\Theta)q(Z) \ln \frac{p(X|\Theta, Z)p(Z | \Theta)p(\Theta)}{q(\Theta)q(Z)} dZ d\Theta \\
&= \int \int q(\Theta)q(Z) \ln p(X|\Theta, Z)  dZ d\Theta  + \int \int q(\Theta)q(Z)\ln p(Z | \Theta) dZ d\Theta \\
&\quad +\int \int q(\Theta) \ln \frac{p(\Theta)}{q(\Theta)}d\Theta +  \int q(Z)\ln\frac{1}{ q(Z)} dZ 
\end{align}

この式において、まずq(\Theta) = {\hat q(\Theta)}という既知の関数として、q(Z)について最大化することを考える。


\displaystyle
\begin{align}
{\mathcal L} \left[{\hat q}(\Theta) q(Z) \right] &=\int \int {\hat q}(\Theta)q(Z) \ln p(X|\Theta, Z)  dZ d\Theta  + \int \int {\hat q}(\Theta)q(Z)\ln p(Z | \Theta) dZ d\Theta \\
&\quad + \int q(Z)\ln\frac{1}{ q(Z)} dZ + {\rm const} \\
&=\int  q(Z) \langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)}  dZ   + \int q(Z)\langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} dZ \\ 
&\quad + \int q(Z)\ln\frac{1}{ q(Z)} dZ + {\rm const} \\
&=\int  q(Z) \ln \exp(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)})  dZ   + \int q(Z)\ln \exp(\langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)}) dZ \\
&\quad + \int q(Z)\ln\frac{1}{ q(Z)} dZ + {\rm const} \\
&= \int q(Z)\ln\frac{\exp \left(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} \right)}{ q(Z)} dZ + {\rm const}
\end{align}

この最後の式をよく見ると、分子をZ積分した時1になれば(つまり確率分布であれば)、これはKLダイバージェンス(に-1をかけたもの)である。そのため、

\displaystyle r(Z) \propto \exp \left(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} \right)
\displaystyle \int r(Z) dZ = 1

と定義すれば、

\displaystyle {\mathcal L} \left[{\hat q}(\Theta) q(Z) \right] = -KL\left[ q(Z) \parallel r(Z) \right] + {\rm const}

になり、下限を最大化するためにはKLダイバージェンスを0にするq(Z)を求めればよい。すなわち、

\displaystyle q(Z) = r(Z) \propto exp \left(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)} \right)

とすればよい。この潜在変数側の更新を変分Eステップという。ちなみに両辺logを取ると、

\displaystyle \ln q(Z) =  \langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(\Theta)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(\Theta)}  + {\rm const}

と書ける。q(Z)がどんな分布かを求める際、このlogを取った形式で計算をすることが多い。

同様の式展開をq(Z) = {\hat q}(Z)という既知の関数としてq(\Theta)について最大化する場合に適用すると、


\displaystyle
\begin{align}
{\mathcal L} \left[q(\Theta) {\hat q}(Z) \right] =
\int q(\Theta)\ln\frac{\exp \left(\langle \ln p(X|\Theta, Z) \rangle_{{\hat q}(Z)} + \langle \ln p(Z | \Theta) \rangle_{ {\hat q}(Z)} \right) p(\Theta)}{ q(\Theta)} d\Theta + {\rm const}
\end{align}

なので、q(\Theta)に関する更新式は、

\displaystyle q(\Theta) \propto \exp \left( \langle \ln p(X | \Theta, Z) \rangle_{{\hat q}(Z)}  +  \langle  \ln p(Z | \Theta) \rangle_{{\hat q}(Z)}  \right) p(\Theta)

となる。こちらは変分Mステップと呼ばれる。こちらも両辺logを取った式を書いておくと、

\displaystyle \ln q(\Theta) =   \langle \ln p(X | \Theta, Z) \rangle_{{\hat q}(Z)}  +  \langle  \ln p(Z | \Theta) \rangle_{{\hat q}(Z)} + \ln  p(\Theta) + {\rm const}

ベイズ推論による教師あり分類 ~ガウス混合分布編~

はじめに

以前、ベイズ推論による機械学習入門という書籍を読み、以下のような記事を公開しました。

yamagensakam.hatenablog.com

この書籍には、多次元ガウス混合モデルによるクラスタリング手法が詳しく述べられていますが、多次元ガウス混合モデルによる教師あり分類の手法は少し触れる程度にとどまっています(クラスタリングで述べられている導出過程を流用すれば、比較的簡単に導出できるためだと思います)。
ということで、今回は多次元ガウス混合モデルによる教師あり分類の導出、実装をやってみました。

ガウス混合モデル

モデル自体は教師ありだろうが、教師なしだろうが同じモデルを考えることができます。まず、全て独立でサンプリングされるD次元のデータが{\bf X}  = \left\{ {\bf x}_1 , \cdots, {\bf x}_N \right\}N個あるとし、そのそれぞれに対応するラベルデータを{\bf S} = \left\{ {\bf s}_1 , \cdots, {\bf s}_N \right\} とします。ここで、{\bf s}1 of K表現と呼ばれ、クラス数がKの場合、所属するクラスに対応する要素のみ1で他は全て0K次元ベクトルです。この{\bf s}をサンプルするための分布として、パラメータ{\bf \pi}を持つカテゴリ分布を考えます。
また、K個のガウス分布が混合したモデルを考えるので、それぞれの分布の平均を{\bf \mu}_1, \cdots, {\bf \mu}_K、精度行列(共分散行列の逆行列)を{\bf \Lambda}_1,\cdots, {\bf \Lambda}_Kとします。なお、各パラメータを個別に扱わない場面では表記を簡単にするため、これらのガウス混合分布のパラメータを全てまとめて{\bf \Theta}と表記します。
更に、カテゴリ分布のパラメータ{\bf \pi}や各ガウス分布{\bf \mu}_k{\bf \Lambda}_kといったパラメータもベイズ推論の枠組みでは、確率分布に従うと考えます。ここでは、それぞれ共役事前分布に従うことにします。具体的には、

\displaystyle p({\bf \pi}) = {\rm dir}({\bf \pi} | {\bf \alpha}) \tag{1}
\displaystyle p({\bf \mu}_k, {\bf \Lambda}_k)  = {\rm NW} ({\bf \mu}_k, {\bf \Lambda}_k | {\bf m}, \beta, \nu, {\bf W} ) \tag{2}

という分布に従うとします。ここで、{\rm dir}(\cdot)はディリクレ分布、{\rm NW}(\cdot)ガウス・ウィシャート分布です(具体的な分布の式は書籍を参照ください)。ここで事前分布として共役事前分布を採用しておくと事後分布や予測分布が解析的に求めることができます。{\bf \alpha}{\bf m}, \beta, \nu, {\bf W} は超パラメータで今回は事前に与えるものとします。また、{\bf m}, \beta, \nu, {\bf W} については、それぞれのガウス分布で異なる値を与えることもできますが、簡単のため今回は全て共通の値とします。
さて、ここまで出てきた確率変数{\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}の同時分布を考えると、以下のようになります。
{
\displaystyle
\begin{eqnarray}
p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}) &=& p({\bf X} | {\bf S}, {\bf \Theta}) p ({\bf S}|{\bf \pi})p({\bf \Theta})p({\bf \pi} ) \\
&=& \left\{\prod_{i = 1}^{N} \prod_{k = 1}^{K} {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1})^{{\bf s}_{i, k}} {\rm Cat}({\bf s}_k | {\bf \pi})\right\} \left\{ \prod_{k = 1}^{K} p({\bf \mu}_k, {\bf \Lambda}_k) \right\} p({\bf \pi}) 
\end{eqnarray}
} \tag{3}

ここで、{\rm N}(\cdot)ガウス分布{\rm Cat}(\cdot)はカテゴリ分布を表し、p({\bf \pi})p({\bf \mu}_k, {\bf \Lambda}_k) は式(1)、式(2)で表される分布です。
このモデルのグラフィカルモデルを下図に示します。この図で{\bf X}{\bf S}の四角で囲まれているのが学習データとして与えられるデータ、{\bf x}_*{\bf s}_*がクラスを予測する対象のデータ(テストデータ)です。


f:id:YamagenSakam:20190210190923p:plain

このように図示しておくと有効分離の考え方を用いて、変数間の条件付き独立性がわかりやすくなり、後々式展開する際に便利です。有効分離については、以下の記事がわかりやすいと思います。

条件付き独立と有向分離を用いた統計モデルの妥当性チェック - StatModeling Memorandum

推論

上述式(3)のモデルから、パラメータの事後分布とそれを元に、新たなデータ{\bf x}_*が与えられた時の{\bf s}_*の確率分布である予測分布を導出します。このモデルの場合、上述のように共役事前分布を事前分布として採用しているため、事後分布および予測分布は解析的に求めることができます。ここでは、逐次的にデータが与えられるとして、オンラインで予測分布を更新していけるような式を導出します。

事後分布の導出

今回扱う問題では{\bf X}{\bf S}がデータとして与えられるので、これらが与えられている下でのパラメータの事後分布を求めます。つまり、

 \displaystyle p({\bf \Theta}, {\bf \pi}|{\bf X}, {\bf S}) = \frac{p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta})}{p({\bf X}, {\bf S})}  \varpropto p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}) \tag{4}

を求めることになります。なお、クラスタリングの場合、{\bf S}がデータとして与えられているわけではないので、求めるべき事後分布としてはp({\bf \Theta}, {\bf \pi}, {\bf S}|{\bf X}) となり、こちらは解析的に求めることが難しくなります(なので、ギブスサンプリングや変分推論などの近似手法を使う)。
ここで式(3)を見ると、\Theta\piは完全に分離していることがわかります。つまり、{\bf X}{\bf S}が与えられた下では、\Theta\piは独立であり、

 \displaystyle  p({\bf \Theta} |{\bf X}, {\bf S}) \varpropto  p({\bf X} | {\bf S}, {\bf \Theta})p({\bf \Theta}) \tag{5}
 \displaystyle  p({\bf \pi} |{\bf X}, {\bf S}) =p({\bf \pi} |{\bf S})  \varpropto   p ({\bf S}|{\bf \pi})p({\bf \pi} ) \tag{6}

を求めればよいことになります。
それでは、式(5)を展開してみます。式(5)の右辺について、対数をとってみると、

{
\displaystyle 
\begin{eqnarray}
\ln p({\bf X} | {\bf S}, {\bf \Theta}) p({\bf \Theta}) &=& \ln \prod_{i = 1}^{N} \prod_{k = 1}^{K} {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1})^{{\bf s}_{i, k}} p({\bf \mu}_k, {\bf \Lambda}_k)\\
 &=& \sum_{k = 1}^{K} \left\{ \sum_{i = 1}^{N} s_{i,k} \ln {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1}) + {\rm NW} ({\bf \mu}_k, {\bf \Lambda}_k | {\bf m}, \beta, \nu, {\bf W} ) \right\}
\end{eqnarray}
} \tag{7}

となります。これを見ると各kに対して分離しているのがわかるので、各kに対し{\bf \mu}_k, {\bf \Lambda}_kそれぞれ独立に事後分布を求めればよいことになります。また、式(7)の中括弧の中はガウス・ウィシャート分布を事前分布として用いたガウス分布パラメータの事後分布推定の計算となります。この計算方法は件の書籍に詳しく載っているので、ここでは割愛しますが、最終的な形は以下のようになります。
\displaystyle p({\bf \mu}_k, {\bf \Lambda}_k | {\bf X}, {\bf S}) =  {\rm WA}({\bf \mu}_k, {\bf \Lambda}_k |  \hat{{\bf m}}_k, {\hat \beta}_k, {\hat \nu}_k, {\hat {\bf W}}_k) \tag{8}
ただし、

{
\displaystyle
\begin{eqnarray}
{\hat \beta}_{k} &=& \sum_{i = 1}^{N} s_{i, k} + \beta  \\
{\hat \nu}_{k} &=& \sum_{i = 1}^{N} s_{i, k} + \nu \\ 
{\hat {\bf m}}_{k} &=& \frac{\sum_{i = 1}^{N} s_{i, k} {\bf x}_{i} + \beta {\bf m}  } {{\hat \beta}_{k}} \\
{\hat {\bf W}}_{k} ^{-1} &=&   \sum_{i = 1}^{N} s_{i, k} {\bf x}_{i} {\bf x}_{i}^T + \beta {\bf m}{\bf m}^T - {\hat \beta}_{k}{\hat {\bf m}}_{k}  {\hat {\bf m}}_{k}^{T} + {\bf W}^{-1}\\
\end{eqnarray}
} \tag{9}

です。同様に式(6)はディリクレ分布を事前分布としてカテゴリ分布の事後分布そのものであるため、

\displaystyle p({\bf \pi} |{\bf S}) = {\rm dir}({\bf \pi}| {\hat {\bf \alpha}}) \tag{10}

というディリクレ分布になります。ただし、

\displaystyle  \alpha_{k} = \sum_{i = 1}^{N} s_{i, k} + \alpha \tag{11}

です。

ここで、式(9)および式(11)の各超パラメータをデータが与えられるたびに逐次更新できる形に変更すると、

{
\displaystyle
\begin{eqnarray}
{\hat \beta}_{n, k} &=&  {\hat \beta}_{n-1, k} + s_{n, k} \\
{\hat \nu}_{n, k} &=&  {\hat \nu}_{n-1, k}  + s_{i, k} \\ 
{\hat {\bf m}}_{n, k} &=&  \frac{s_{n, k} {\bf x}_{n} + {\hat \beta}_{n-1, k} {\hat {\bf m}}_{n-1, k}  } {{\hat \beta}_{n, k}}\\
{\hat {\bf W}}_{n, k}^{-1} &=& { \hat {\bf W}}_{n-1, k}^{-1} + s_{i, k}{\bf x}_{i} {\bf x}_{i}^T +{\hat \beta}_{n-1, k}{\hat {\bf m}}_{n-1, k}  {\hat {\bf m}}_{n-1, k}^{T}  - {\hat \beta}_{n, k}{\hat {\bf m}}_{n, k}  {\hat {\bf m}}_{n, k}^{T} \\
{\hat \alpha}_{n, k} &=&  {\hat \alpha}_{n-1, k}  + s_{i, k} 
\end{eqnarray}
} \tag{12}

となり、これによりn-1の時点で、計算した超パラメータを用いて、n番目のデータが来たときの超パラメータに更新することが可能です(初期値は事前分布の超パラメータ{\bf m}, \beta, \nu, {\bf W}, \alpha )。

予測分布の導出

最終的に欲しいのは、新たなデータ{\bf x}_*が与えられた時のクラス分類結果{\bf s}_*です。そのため、{\bf x}_*が与えられた下での{\bf s}_*の予測分布を求めます。予測分布は、

\displaystyle p({\bf s}_* | {\bf x}_*, {\bf X}, {\bf S}  ) = \frac{p({\bf s}_* , {\bf x}_*, {\bf X}, {\bf S}  )}{p( {\bf x}_*, {\bf X}, {\bf S} )} = \frac{p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S}  )p( {\bf X}|{\bf s}_*, {\bf S})p({\bf s}_*|{\bf S})p({\bf S})}{{p( {\bf x}_*, {\bf X}, {\bf S} )}}  \tag{13}

となります。ここで、{\bf s}_*に関係のない項は無視していいので、p({\bf S})p( {\bf x}_*, {\bf X}, {\bf S} )は考慮から外します。更に、p( {\bf X}|{\bf s}_*, {\bf S})について、グラフィカルモデルから有効分離性を読み解くと{\bf S}が与えられた下では、{\bf X}{\bf s}_*は独立であり({\bf X}に繋がる{\bf S}が観測され、かつ、head-to-head型を形成している{\bf s}_*, {\bf x}_*, {\bf \Theta} において{\bf x}_*が観測されていないため)、p( {\bf X}|{\bf s}_*, {\bf S})= p( {\bf X}| {\bf S}) となるため、こちらも考慮しなくてよくなります。まとめると、式(13)は

\displaystyle p({\bf s}_* | {\bf x}_*, {\bf X}, {\bf S}  ) \varpropto p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S}  ) p({\bf s}_*|{\bf S}) \tag{14}

と、2つの項のみ考えればよくなります。1つ目の項は、
 \displaystyle p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S}  ) = \int  p({\bf x}_* | {\bf s}_* ,{\bf X},{\bf S}, {\bf  \Theta} )  p({\bf  \Theta}| {\bf s}_*,{\bf X},{\bf S}) d {\bf \Theta} \tag{15}

であり、ここでもグラフィカルモデルから有効分離性を考えるとp({\bf x}_* | {\bf s}_* ,{\bf X},{\bf S}, {\bf  \Theta} )  =p({\bf x}_* | {\bf s}_* , {\bf  \Theta} ) p({\bf  \Theta}| {\bf s}_*,{\bf X},{\bf S})  = p({\bf  \Theta}| {\bf X},{\bf S})とすることができます。従って、

\displaystyle p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S}  ) =  \int  p({\bf x}_* | {\bf s}_* , {\bf  \Theta} )  p({\bf  \Theta}| {\bf X},{\bf S}) d {\bf \Theta} \tag{16}

となります。
 p({\bf  \Theta}| {\bf X},{\bf S}) は、式(5)、式(8)、式(9)に記している事後分布そのものです。今、各kに対して{s}_{*,k} = 1であるそれぞれの予測分布を求めたいので、式(16)より

\displaystyle p({\bf x}_*| {\bf s}_{*,k} = 1 , {\bf X}, {\bf S}  ) =  \int  \int p({\bf x}_* | {\bf  \mu}_k, {\bf \Lambda}_k )  p({\bf  \mu}_k, {\bf \Lambda}_k | {\bf X},{\bf S}) d {\bf  \mu}_k d {\bf  \Lambda}_k \tag{17}

を求めればよいことになります。これは、ガウス・ウィッシャート分布を事前分布として用いた場合のガウス分布の予測分布計算になります。こちらも予測分布導出は書籍に載っているので、導出過程は省略し最終形を示します。

\displaystyle p({\bf x}_*| {\bf s}_{*,k} = 1 , {\bf X}, {\bf S}  ) = {\rm St}({\bf x}_* |{\hat {\bf m}}_{n, k}, \frac{(1-D+ {\hat \nu}_{n, k}){\hat \beta}_{n, k} }{1 + {\hat \beta}_{n, k} } {\hat {\bf W}}_{n,k}, 1-D+ {\hat \nu}_{n, k}) \tag{18}

ただし、 {\rm St}(\cdot)は多次元のスチューデントのt分布で、D{\bf x}の次元です。
次に、式(14)のp({\bf s}_*|{\bf S})ですが、こちらも同様に、

\displaystyle p({\bf s}_*|{\bf S}) = \int p({\bf s}_*|{\bf S},{\bf \pi})p({\bf \pi}| {\bf S}) d {\bf \pi} =  \int p({\bf s}_*|{\bf \pi})p({\bf \pi}| {\bf S}) d {\bf \pi} \tag{19}

となります。ただし、第2式から第3式への展開はグラフィカルモデルより {\bf \pi}を観測した下での{\bf S} {\bf s}_*の条件付き独立性を利用しました。そしてこの式(19)は、ディリクレ分布を事前分布として使ったカテゴリ分布の予測分布です。ということで、この予測分布についても導出過程は割愛し最終形を示します。

\displaystyle p({\bf s}_*|{\bf S})  = {\rm Cat}({\bf s}_*|{\hat {\bf \eta}}_{n}) \tag{20}

ただし、{\hat {\bf \eta}_{n}}の各要素{\hat \eta}_{n,k}
\displaystyle {\hat \eta}_{n,k}=  \frac{{\hat \alpha}_{n, k}}{\sum_{i=1}^K  {\hat \alpha}_{n, i}}  \tag{21}
です。
ここまでの話をまとめると、n番目のデータが与えら手た時の予測分布の更新式は、式(14)、式(17)、式(20)、式(21)より

 \displaystyle p(s_{*, k} = 1 | {\bf x}_*, {\bf X}, {\bf S}  )  \varpropto  {\hat \eta}_{n,k} {\rm St}({\bf x}_* |{\hat {\bf m}}_{n, k}, \frac{(1-D+ {\hat \nu}_{n, k}){\hat \beta}_{n, k} }{1 + {\hat \beta}_{n, k} } {\hat {\bf W}}_{n,k}, 1-D+ {\hat \nu}_{n, k}) \tag{22}

となります(各超パラメータは式(12)、式(21))。そして、各クラスの所属確率の平均を得たいときは、k = 1 \cdots Kに対し、式(22)の右辺に{\bf x}_*を代入し、得られる値を合計が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分布のクラスで、この実装については、以下の記事の内容を参考にさせてもらいました。

<ベイズ推論> 多次元ガウス分布の学習 - Helve’s Python memo

実験・結果

ということで、以下のようなコードを組んで実験してみました。クラス数はK=3で、学習データ数はそれぞれN_1=50, N_2 = 60, N_3 = 70と少し偏りを持たせています。
また、せっかく事後分布を逐次学習できるようにしたので、学習データが与えられるたびに予測平均がどのように変化するかをアニメーションにしてみました。

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()

そして、こちらが結果です。

f:id:YamagenSakam:20190210212509g:plain

この図で、各色がそれぞれのクラスを表し、追加されていく各点が学習データを表しています。また、背景色が予測分布が算出した各座標点の各クラスへの所属確率を表しています。
特に左上に辺りの領域に着目すると、最初の方は青色のクラスと予測されていましたが、青と緑のクラスのデータを観測していくにつれ、緑色のクラスだと予測されるようになっていくのがわかります。これは青色のクラスが左下から右上にかけての共分散が大きいのに対し、緑色のクラスの共分散が左上から右下への共分散が大きく、その特徴をサンプルが増えるにつれ、正しく推論できるようになっていくためだと考えられます。

まとめ

今回、ガウス混合分布を用いた教師あり分類をベイズ推論の枠組みを適用して実装してみました。教師ありの分類の場合は、ガウス混合分布モデルによる予測は解析的に求めることができます。
ベイズ推論の枠組みで分類を行う場合は他にもロジスティック回帰をベイズの枠組みで定式化する方法などがあります。今後は、このベイズロジスティック回帰を含めもう少し複雑なモデルを学び実装していきたいと思います。

ベイズ推論による機械学習入門を読んだので軽くメモ

はじめに

タイトルの通りです。以下の本を一通り読んだので、再度読む際の手助けになるようメモを残します。

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

まだ理解が怪しい部分も多々ありますが、取りあえず忘れる前に。。。

全体を通してのポイント

ベイズ推論の枠組み

ベイズ推論においては、観測データも分布のパラメータ、回帰等の重み係数も潜在変数も確率変数として扱う。観測されたデータをDとし、未知の確率変数(分布のパラメータ、重み係数、潜在変数等々)を{\bf \theta}とした時、以下の2ステップを行うのがベイズ推論の枠組み。

  1. モデル\displaystyle p(D, {\bf \theta})の構築(モデリング
  2. 事後分布\displaystyle p( {\bf \theta} | D) = \frac{p(D, {\bf \theta})}{p({D})}の計算

この事後分布を求めることがベイズ推論における「学習」。ベイズ推論は全てこの枠組みの中で扱われるので、後は如何にモデリングするかと、如何に求めたい確率変数の事後分布を計算するか、という話になる。モデルが複雑な場合、事後分布は解析的に求められない場合もあるので、その時は変分推論やギブスサンプリングなどを使って近似的に事後分布を求める。

予測分布

更に、求めた事後分布を使って、未観測のデータ{\bf x}_*に対し、何らかの知見を得たい場合もある(知見の内容はタスクによって異なる。例えば、回帰の場合は、入力{\bf x}_*に対応する出力y_*を求める等)。その時は以下の予測分布を計算する。


\displaystyle p({\bf x}_* | D) = \int p( {\bf x}_* | {\bf \theta}) p( {\bf \theta} | D) d  {\bf \theta}

1章

基本的にこの章は後の章で語られる内容のアウトラインが載っている。この章で語られる内容で特に重要なのは、上記のベイズ推論のポイントと以下の関係式。

  • 周辺化 \displaystyle p(y) = \int p(x, y) dy
  • 条件付き分布 \displaystyle p(x | y) = \frac{p(x, y)}{p(y)}

また、独立の諸概念も重要。

  • 変数間の独立性p(x, y) = p(x)p(y)(これが成り立つ時、xyは独立という)
  • 条件付き独立p(x,z|y) = p(x|y)p(z|y)yを観測した下でxzは独立だが、yを観測しなければxzは独立でないことに注意)

2章

確率密度関数については書籍を参照するとして、取りあえず期待値が大事。エントロピーとかKLダイバージェンスとかも大事だけど、期待値がわからないと後々お話にならなくなる。特に、


\displaystyle \langle  f({\bf x}) \rangle_{p({\bf y})} = \int f({\bf x}) p({\bf y}) d{\bf y} =  f({\bf x})  \int p({\bf y}) d{\bf y}  =  f({\bf x })
\displaystyle \langle  f({\bf x}) \rangle_{p({\bf x}, {\bf y})} = \int \int f({\bf x}) p({\bf x}, {\bf y}) d{\bf x}d{\bf y} =  \int f({\bf x}) \int  p({\bf x}, {\bf y}) d{\bf y} d{\bf x}  = \langle f({\bf x }) \rangle_{p({\bf x})}

と、期待値を取る関数に出てこない確率変数は、期待値計算から外せることは重要。
また、同時分布の期待値とか、条件付き期待値とかは込み入ってくると混乱するので注意。迷った時は期待値の定義、


\displaystyle \langle   f({\bf x}) \rangle_{p({\bf x})}  = \int f({\bf x}) p({\bf x}) d{\bf x}

に立ち返るのが吉。

確率密度関数の対数を取った関数\ln p({\bf x})、基本的な期待値( \langle  {\bf x} \rangleとか \langle \ln {\bf x} \rangleとか)も、後々の章(特に4章、5章の変分推論辺り)で何度も使う。本当に何度も使うので、何度も見返す。

3章

2章で挙がっている確率分布の各パラメータおよび線形回帰の重み係数をベイズ推論の枠組みで推定する章。事前分布として共役事前分布を使うので、(この章で示されている問題については)事後分布が解析的に計算でき、しかも事前分布として同じ形の分布で求まるので、データを得るたびに逐次的に更新することも可能。考え方としては、

  1. 学習したいパラメータ\thetaの事前分布p(\theta)として共役事前分布を設定
  2. こうすると事後分布は事前分布と同じ形になり、p( {\bf \theta} | D) = \frac{p(D, {\bf \theta})}{p({D})}= \frac{p(D|{\theta})p(\theta)}{p({D})} \propto p(D|{\theta})p(\theta) としてp(D)は定数として扱うことが可能(後から帳尻を合わせる)。
  3. 尤度p(D | \theta)p(D | \theta) = \prod_{i=1}^N p(x_i|\theta)と観測データが独立だと仮定(後の章で観測データ間にも独立ではないモデルもあり)

また、計算のテクニックとして、以下のようなものが使われている。

  • 事後分布が事前分布と同じ形になることから、\ln p(D|{\theta})p(\theta)をうまく展開し、2章で計算済みの事後分布の対数をとった式のパラメータと対応を取ることで、事後分布のパラメータを計算
  • 予測分布を求める際、積分の計算が難しい時は、ベイズの定理に立ち返りp({\bf x}_*)=\frac{p({\bf x}_* | \theta) p(\theta)}{p(\theta | {\bf x}_* )}となることから、上記の事後分布計算と同じく対数をとって分布の式とパラメータの対応を取ればよい
  • 予測分布計算時、いきなり事後分布のパラメータで予測分布求めようとすると、ややこしくなるケースもあり、その場合はまず事前分布を用いて予測分布を計算し、後から対応するパラメータを事後分布のものに置き換えればよい

4章

生成モデル

観測データの生成過程を仮定して、確率分布p(D, {\bf  \theta})モデリングしたものが生成モデル。例えば、観測されたデータに多峰性が見れられる場合、これを単一の分布のみで表現しようとしても、うまく表現しきれない。このような多峰性のあるデータは、

  1. K個のクラスタがあり、観測データ{\bf x}はそのいずれかのクラスタ(例えば、\theta_kを持つk番目のクラスタ)から生成されたと仮定p(\bf{x} | \theta_k)
  2. {\bf x}が割り当たっているクラスタ{\bf s}  \in \left\{0, 1 \right\}^K \ \ {\rm s.t.} \ \ \sum_{i=1}^K s_i = 1(割り当たっているクラスタの要素のみ1、それ以外0)で表現すると、{\bf x}の確率分布はp({\bf x}| {\bf s} ,\Theta) = \prod_{k=1}^K p(\bf{x} | \theta_k)^{s_k}
  3. このクラスタの割り当て{\bf s}も割り当て比率パラメータ\piに支配されるカテゴリ分布の確率変数と仮定p({\bf s}| \pi)
  4. クラスタのパラメータの事前分布p(\Theta)を仮定(ここでも共役事前分布使うと幸せになれる)
  5. クラスタ割り当てのカテゴリ分布の事前分布p(\pi)を仮定(こちらもカテゴリ分布の共役事前分布であるディリクレ分布を使うと幸せ)
  6. 以上のことから、データ生成のモデルをp({\bf X}, {\bf S}, \Theta, \pi) = p({\bf X}|{\bf S} , \Theta) p({\bf S}| \pi)  p(\Theta) p(\pi) = \prod_{i = 1}^N p({\bf x}_i| {\bf s}_i ,\Theta) p({\bf s}_i| \pi)  p(\Theta)p(\pi)と仮定

等のようにモデリングする。
このモデルから、観測データに対するクラスタ割り当ての事後分布\displaystyle p({\bf S}|{\bf X}) = \int \int \frac{p({\bf X}, {\bf S}, \Theta, \pi) }{p({\bf X})} d\Theta d\pi を計算したいが、モデルが複雑すぎてp({\bf X})が解析的に求まらない。そのため、近似的に事後分布の推定を行う。

ギブスサンプリング

上のモデルを例にすると、

  1. \Theta, \piが既に観測されているものとして(2でサンプリングした\Theta, \pi を使う)、{\bf S} \sim p({\bf S}| {\bf X}, \Theta, \pi)とサンプリング
  2. {\bf S}が既に観測されているものとして(1でサンプリングした{\bf S}を使う)、\Theta, \pi \sim p(\Theta, \pi | {\bf X}, {\bf S} )とサンプリング
  3. この操作を繰り返す

この繰り返し回数が十分に多い場合、真の分布に近づくことが知られているため、繰り返し続けると所望のデータをサンプリングできるようになる。もちろん、サンプリングする対象となる分布p({\bf S}| {\bf X}, \Theta, \pi)p(\Theta, \pi | {\bf X}, {\bf S} )は十分にサンプリングしやすい分布である必要があるので、式を展開して2章で挙がっているような基本的な分布を導かないといけない(モデルによってはそれができないこともあるので、共役事前分布を使うと幸せになるというのはここに繋がる)。

変分推論


\displaystyle p({\bf S}, \Theta, \pi | {\bf X}) \simeq q({\bf S}) q(\Theta, \pi)

等のように、各確率変数間の独立性を仮定した分布で近似するのが変分推論。真の分布とのKLダイバージェンスが最小となるq({\bf S}) q(\Theta, \pi)を求めればよく、それは


\displaystyle \ln q({\bf S})= \langle \ln p({\bf X}, {\bf S}, \Theta, \pi ) \rangle_{q(\Theta, \pi)}
\displaystyle \ln q(\Theta, \pi) = \langle \ln p({\bf X}, {\bf S}, \Theta, \pi ) \rangle_{q({\bf S})}

として計算できる(あくまで、上記モデルに対する近似式。もう少し一般的な表現は書籍参照)。こちらも一方の結果を用いで、もう一方の計算を行う形になっているため、更新を十分な回数繰り返す必要がある。また、実際には、下記のような式展開して期待値を解析的に求める。

  1. 変分推論を基に、各近似分布q(\cdot)の形を求める。2章に述べられている期待値の性質(上記の無関係な確率変数が期待値計算から外せる性質含む)をうまく使い、q(\cdot)が何分布の関数かがわかるまで式を展開する。
  2. 全ての近似分布の形がわかったら、それぞれの期待値を計算する。これまた2章で計算済みの各分布の基本的な期待値を流用する。

ここでもギブスサンプリングと同様、期待値を解析的に計算するために、2章で挙がっているような基本的な分布を導く必要があるので、やはり共役事前分布を使った方がハッピー。

崩壊型ギブスサンプリング

上記のギブスサンプリングでは、パラメータ\Theta, \pi についてもサンプリングで求めたが、


\displaystyle p({\bf X}, {\bf S}) = \int \int p({\bf X}, {\bf S}, \Theta, \pi) d\Theta d\pi

として周辺化すれば、p({\bf S}| {\bf X})からサンプリングするだけでよくなる。ただし、周辺化した結果が十分にサンプリングしやすい確率分布になることが条件。このモデルの例では、周辺化した結果{\bf s}_1 \cdots, {\bf s}_Nが互いに依存し合うようになり(元々は{\bf X},\Theta, \piを与えられた下での条件付き独立)、このままではN個の潜在変数を同時にサンプリングしなければならず計算量が莫大になる。この問題をp({\bf s}_n | {\bf X}, {\bf S}_{ \backslash n})と、サンプリングしたい{\bf s}_n以外の潜在変数{\bf S}_{ \backslash n}は全て観測済みとして、個別にサンプリングするアイデアで解決。

5章

各種応用に対して、モデルの構築と変分推論(一部ギブスサンプリング)を導出している。この中のいくつかは実際に実装してみようと思うので、ここでは考え方のみメモ。

線形次元削減

観測されたデータ{\bf y}が低次元の潜在変数{\bf x}を元に、 {\bf W}^T {\bf x} + {\bf \mu}という線形変換で生成されると考える。また、実際にはノイズが乗っているはずなので、ガウス分布N({\bf y}| {\bf W}^T {\bf x} + {\bf \mu}, \sigma^2 {\bf I})で生成されると仮定。更に他のパラメータについても事前分布がガウス分布に従うと仮定(\sigma^2は既知のものとして話を進めているが、ガンマ事前分布など使って推定してもよい)。後は変分推論の期待値計算をすれば低次元の潜在変数{\bf x}や変換行列 {\bf W}を推測することができる。

非負値行列因子分解

観測された行列{\bf X}が全ての要素非負であるとし、{\bf X} \simeq {\bf W} {\bf H}という非負値の行列に近似分解することを考える。ここで、S_{dmn} \simeq W_{dm}H_{mn}となる潜在変数{\bf S}を考えることで、効率的な変分推論が導き出せる。また、非負値なので{\bf S}ポアソン分布、{\bf W}, {\bf H}の事前分布はガンマ分布を考える。

隠れマルコフモデル

これまで観測データ間はたがいに独立して得られるものとして扱っていたが、時系列データ等ではその仮定は成り立たない。隠れマルコフモデルでは、

  1. 観測されるデータ{\bf x}_nは現在の状態{\bf s}_nに基づいてデータが観測される(状態{\bf s}_nクラスタ割り当ての潜在変数として持つ混合分布)。
  2. データを観測したら、遷移確率に基づいて次の状態に遷移する({\bf s}_n \rightarrow {\bf s}_{n + 1})。

というモデルを構築する。変分推論は状態の潜在変数{\bf S}の分布を\prod_i^N q({\bf s}_i)と分離して近似する完全分解変分推論とq({\bf S})と分離せずに近似する構造化変分推論がある(その中間であるミニバッチを使った近似もあり)。前者は特別なことをせずに、変分推論の期待値を導き出せるが、後者はうまいこと式変形してフォワードバックワードアルゴリズムと呼ばれる動的計画法が使える形に持っていかないと、計算量が増大する。

トピックモデル

複数の文書を観測データとして、単語の出現頻度を元に文書を解析。各文書中には複数の潜在トピックが存在すると仮定。文書d中のn番目に出てくる単語{\bf w}_{dn}には、どのトピックから発生したかというトピック割り当て潜在変数{\bf z}_{dn}が紐づき、各単語は対応するトピックkのカテゴリ分布に従うと仮定。また、トピック割り当て自体もカテゴリ分布に従うと仮定し、これらのカテゴリ分布のパラメータはディリクレ事前分布に従うモデルを構築。後は頑張って変分推論の期待値計算。

テンソル分解

協調フィルタリングに時間軸を導入し、その時のトレンドを踏まえた推薦を行うようにしたモデル。時間方向の軸には前回の状態を平均とするガウス分布を仮定(線形動的システム)。ここでは、完全分解変分推論を考えているので、推論計算自体はこれまでと同様。

ロジスティック回帰、ニューラルネットワーク

ロジスティック回帰、ニューラルネットワークベイズ推論版。これまでの近似(平均場近似)とは違い、ガウス分布で事後分布を近似。変分パラメータを用いた勾配法で真の分布とのKLダイバージェンス最小化を行う。この勾配法を使うため、再パラメータ化トリックを使うが、これが良く分からない。

おわりに

兎にも角にもわかりやすいです。式展開がわからなくなって何度も行ったり来たりしましたが、PRML10章で卒倒した僕でも割りと短期間で最後まで行けたので、ベイズ推論の入り口として、適切な書籍の1つだと思います。ただ機械学習自体の入門としては少ししんどいので、機械学習全体を俯瞰したのち、その1つの枠組みとしてベイズ推論があると捉えられるようになった上で、読むべきかなと思います。

2019年の目標

はじめに

あけましておめでとうございます。今年は『有言実行』というものに挑戦してみようと思い、まず今年の目標を立てそれを本記事に晒します*1。後半雑になっていきますが、目標とそれを達成するためのアプローチを書いたつもりです。

競プロ

  • atcoder青になる
    • 動的計画法に弱いので典型問題を解く
    • 300点問題を解くのが遅いので、速度重視で300穴埋め
    • 400~600穴埋め
    • 主要言語をpythonからC++に乗り換え
    • セグメント木、ネットワークフロー系を理解

数学・機械学習

その他

  • 部屋をきれいに保つ
    • 週1で掃除機をかける
    • 定期的にゴミ出しをする(特にペットボトルを溜め込まない)
    • ゴミ袋をケチらない
  • 体重5kg減
    • 1日の歩数10000歩以上(現状維持)
    • 週1でジム(ランニング・筋トレ中心から水泳中心に変更?)
    • カロリーを意識(朝:400kcal以内、昼:600kcal以内、夜:800kcal以内)
    • 春からジョギング再開
  • DTM始める
    • とりあえずやってみる
  • ブログの書きかけ記事消化(タイトルだけ書いてるものを合わせると15記事ほど)
    • 気が向いたやつから
  • 行き詰ったら、すぐtwitter、slackを見る癖を直す
    • どうすれば・・・

興味あるけどやれてなくて、たぶん今年もやらないだろうなぁと思うやつ

  • ゲームAIコンテスト
    • たぶんこの中だと1番やる可能性高い
  • ラソン(競プロ)
    • これもやる可能性が高い方
  • ラソン(物理)
    • 数年前からハーフに出たいなーとか思いつつ、年々体力が減少して今や10kmすら完走が危うい・・・
  • 量子情報
    • 社内で勉強会とかがあればやるかも
  • 現代制御
  • 離散最適化、離散凸解析
  • 強化学習
  • スマホアプリ開発
  • 婚活
  • 楽器(ベース辺り)
    • ベースは数年前に衝動買いしたから機材はある
  • 格闘技
    • 運動音痴おじさんでも始められるヤツがあれば・・・
  • 資産運用
    • 証券口座は開設したからなんかやるかー

そういえば・・・

確かこのブログでは言っていなかったと思いますが、僕も昨年8月に電機メーカーを退職し、9月よりIT企業にて速いアルゴリズム考えたり、速いプログラム書いたりするお仕事に従事しています。

*1:仕事の目標は10月に職場で立てたので、プライベートの目標のみ

競プロとかに使うアルゴリズム実装メモ(二分探索、2次元累積和、しゃくとり法)

はじめに

アルゴリズムメモ第3段です。今回は二分探索法、2次元累積和、しゃくとり法と様々な問題に使える汎用的なアルゴリズムを書いていきます。
今回は勉強のため、アルゴリズムの本質的な部分を記述した抽象クラスと実際の問題を解く具象クラス(関数オブジェクト)を分け、テンプレートメソッドパターンもしくはストラテジパターンを用いてコーディングしました。

二分探索

用途

  • 広義単調増加関数または広義単調減少関数f(x)が条件を満たし、かつ、条件成立と不成立の境界となる点xを探索(条件成立と不成立の境界が1つのみ)
  • ソート済み配列から条件が成立する値が存在するかの判定

アルゴリズムのポイント

  1. 変数lを探索範囲の最小値、変数rを最大値で初期化
  2. mr,lの中点(例えば、m = (l + r)/2で求める)として、f(m)が条件を満たすか否かを判定
    • このとき、見つけたい点(条件成立と不成立の境界)がmより小さいか、大きいかはf(x)が広義単調増加(減少)関数なのでわかる
  3. 見つけたい点がm以上ならl = m、そうでなければr = mとする
  4. この2~3の操作をr - l \leq  \epsilon  を満たすまで繰り返す(\epsilonは収束判定基準の定数で、問題に応じて適切な値を定める)
  5. 収束したら、その時点のl,rから最終的な出力値を決める

抽象クラスの実装

上記アルゴリズムを実現する抽象クラスを作成し、見つけたい点が現時点以上かの判定、条件判定、中点計算、戻り値計算を純粋仮想関数とし、上記アルゴリズムをテンプレートメソッドパターンで実装。

#include <vector>

template <typename T>
class abstract_binary_search
{
private:
    virtual bool is_in_right(T x) = 0; //求めているのがxより大きいか
    virtual bool is_converge(T L, T R) = 0; //収束判定
    virtual T ret_val(T L, T R) = 0; //見つからなかった場合に返す値
    virtual T get_middle(T L, T R) = 0; //中央値を返す

protected:
    T search_rec(T L, T R);

public:
    virtual T search() = 0;
    virtual bool is_match(T x) = 0;  //それが求めている値か

    virtual ~abstract_binary_search();
};

template <typename T>
T abstract_binary_search<T>::search_rec(T L, T R) {
    if (is_converge(L, R) == true) {
        return ret_val(L, R);
    }

    T mid = get_middle(L, R);

    if (is_in_right(mid) == true) {
        return search_rec(mid, R);
    }

    return search_rec(L, mid);

}


template <typename T>
abstract_binary_search<T>::~abstract_binary_search()
{
}

具体的な問題を解く実装

Lower bownd

素数nの数列が与えられて、指定された値q以上の値が入っているindexを求める問題。もし数列内にq以上の値がなければnを返す。

#include <iostream>

class lower_bound :
    public abstract_binary_search<int>
{
private:
    std::vector<int> data_list;//検索する数列
    int N; //数列のサイズ
    int search_data;//検索対象データ
    virtual bool is_in_right(int x); //求めているのがxより大きいか
    virtual int ret_val(int L, int R); //見つからなかった場合に返す値
    virtual bool is_converge(int L, int R); //収束判定
    virtual int get_middle(int L, int R);
public:
    lower_bound(const std::vector<int>& data_list);
    void set_search_data(int a_search_data);
    virtual int search();
    virtual bool is_match(int x);  //それが求めている値か
    ~lower_bound();
};

bool lower_bound::is_in_right(int x) {
    return data_list[x] <= search_data;
}


int lower_bound::ret_val(int L, int R) {
    if (search_data <= data_list[L])
        return L;
    else
        return L + 1;
}

int lower_bound::get_middle(int L, int R) {
    return (L + R) / 2;
}

lower_bound::lower_bound(const std::vector<int>& data_list):
    data_list(data_list)
{
    N = data_list.size();
}

int lower_bound::search() {
    return search_rec(0, N);
}

bool lower_bound::is_match(int x) {
    if (x >= N) return false;
    return data_list[x] == search_data;
}

bool lower_bound::is_converge(int L, int R) {
    return R - L <= 1;
}

void lower_bound::set_search_data(int a_search_data) {
    search_data = a_search_data;
}

lower_bound::~lower_bound()
{
}

int main()
{

    int n, q;
    std::cin >> n;
    std::vector<int> S(n);
    
    for (int i = 0; i < n; i++) {
        int  s;
        std::cin >> s;
        S[i] = s;
    }
    std::cin >> q;

    int count = 0;
    
    lower_bound lbs(S);
    for (int i = 0; i < q; i++) {
        int  t;
        std::cin >> t;
        lbs.set_search_data(t);
        
        int res = lbs.search();
        std::cout << res << std::endl;
        if (lbs.is_match(res)) 
            count++;
    }
    std::cout << count << std :: endl;
    return 0;
}
POJ 1064:Cable master

蟻本に二分探索の例として載っている問題。問題はリンクを参照。xが切り取る紐の長さ、f(x)が作れる同じ長さ(長さ=x)の数だとすると、f(x)は広義単調減少関数なので二分探索法が使える。
つまり、f(x)が条件を満たしつつ、xが最大となる点を二分探索する(条件を満たさなければxを小さくし、満たすならまだ長くできると考えxを大きくする)。

#include <vector>
#include <algorithm>
#include <functional>
#include <iostream>
#include <iomanip>
#include <math.h>

class cable_master :
    public abstract_binary_search<double>
{
private:
    double tol;
    std::vector<double> len_list;
    int K;
    int N;
    bool is_condition(double x);
    virtual bool is_in_right(double x); //求めているのがxより大きいか
    //virtual bool is_in_left(double x); //求めているのがxより小さいか
    virtual double ret_val(double L, double R); //見つからなかった場合に返す値
    virtual double get_middle(double L, double R);
    virtual bool is_converge(double L, double R); //収束判定
public:
    cable_master(const std::vector<double>& len_list, int K, double tol);
    virtual double search();
    virtual bool is_match(double x);
    ~cable_master();
};

bool cable_master::is_condition(double x) {
    int count = 0;
    for (int i = 0; i < N; i++) 
        count += int(len_list[i] / x);

    if (count >= K)
        return true;

    return false;
}

bool cable_master::is_in_right(double x) {
    return is_condition(x);
}

double cable_master::ret_val(double L, double R) {
    return R;
}

double cable_master::get_middle(double L, double R) {
    return (L + R) / 2.0; 
}

bool cable_master::is_converge(double L, double R) {
    return R - L < tol;
}

cable_master::cable_master(const std::vector<double>& len_list, int K, double tol) :
    len_list(len_list),
    K(K),
    tol(tol)
{
    N = len_list.size();
}

double cable_master::search(){
    return search_rec(0.0, *std::max_element(len_list.begin(), len_list.end()));

}

bool cable_master::is_match(double x) {
    return false;
}

cable_master::~cable_master()
{
}



int main()
{

    std::vector<double> cable_list;
    int N, K;
    std::cin >> N >> K;

    for (int i = 0; i < N; i++) {
        double len;
        std::cin >> len;
        cable_list.push_back(len);
    }
    cable_master cm(cable_list, K, 0.0001);
    std::cout << std::fixed << std::setprecision(2) << floor(cm.search() * 100) / 100.0 << std::endl;
}

2次元累積和

用途

  • \displaystyle \sum_{i = p_x}^{q_x} \sum_{j = p_y}^{q_y}f(x_i, y_j)の計算

アルゴリズムのポイント

  1. あらかじめ、全てのn,mについて\displaystyle  S_{nm} = \sum_{i = 1}^{n} \sum_{j = 1}^{m} f(x_i, y_j)  を求めておく
    • これは、\displaystyle S_{nm} =f(x_n, y_m) + S_{n-1, m} + S_{n, m-1} - S_{n-1, m - 1}で効率的に計算できる
  2. 目的の値は\displaystyle  \sum_{i = p_x}^{q_x} \sum_{j = p_y}^{q_y}f(x_i, y_j) =  S_{q_x, q_y} - S_{q_x-1, q_y} - S_{q_x, q_y-1} + S_{q_x - 1, q_y-1}  で求まる

抽象クラスの実装

このアルゴリズムはストラテジーパターンで実装する。具体的には、f(x_i, y_j)を計算する関数オブジェクト受け取り、それに応じた部分和を求める。

#include <vector>
#include <functional>

template <typename T>
class func_sum
{
private:
    std::vector<std::vector<T>> sum_table;
protected:
    virtual T acc_sum(int i, int j);
public:
    func_sum(std::function<T(int, int)> func, int nx, int ny) ;
    ~func_sum();
    //目的の区間挿話を求める関数
    T calc_section_sum(int px, int qx, int py, int qy);
};

template <typename T>
T func_sum<T>::acc_sum(int i, int j) {
    return sum_table[i][j];
};

template <typename T>
func_sum<T>::~func_sum() {}

template <typename T>
func_sum<T>::func_sum(std::function<T(int, int)> func, int nx, int ny) {
    sum_table = std::vector<std::vector<T>>(nx + 1, std::vector<T>(ny + 1, T()));
    //j=0の場合の総和の計算
    for (int i = 1; i <= nx; i++) {
        sum_table[i][0] += func(i - 1, 0);
    }
    //i=0の場合の総和の計算
    for (int j = 1; j <= ny; j++) {
        sum_table[0][j] += func(0, j - 1);
    }
    //総和の計算
    for (int i = 1; i <= nx; i++) {
        for (int j = 1; j <= ny; j++) {
            sum_table[i][j] += func(i - 1, j - 1) +
                sum_table[i - 1][j] +
                sum_table[i][j - 1] -
                sum_table[i - 1][j - 1];

        }
    }
}

template <typename T>
T func_sum<T>::calc_section_sum(int px, int qx, int py, int qy) {
    return acc_sum(qx, qy) - 
           acc_sum(qx, py - 1) - 
           acc_sum(px - 1, qy) + 
           acc_sum(px - 1, py - 1);
}

具体的な問題を解く実装

\sum_{i = p_x}^{q_x} \sum_{j = p_y}^{q_y} \frac{x_i}{y_i}の計算

f(x_i, y_j)= \frac{x_i}{y_i}  として、上記アルゴリズムに関数オブジェクトを渡す。

int main()
{
    int N, M;
    int px, py, qx, qy;
    std::cin >> N >> M ;
    std::cin >> px >> qx >> py >> qy;
    std::vector<double> x_list(N, 0);
    std::vector<double> y_list(M, 0);
    for (int i = 0; i < N; i++) {
        double x;
        std::cin >> x;
        x_list[i] = x;
    }
    for (int i = 0; i < M; i++) {
        double y;
        std::cin >> y;
        y_list[i] = y;
    }
    auto func = [](int i, int j, 
                   const std::vector<double>& x_list,
                   const std::vector<double>& y_list) { return x_list[i] / y_list[j]; };

    std::function<double(int, int)>  func_curried = std::bind(func, 
                                                              std::placeholders::_1, 
                                                              std::placeholders::_2, 
                                                              x_list, 
                                                              y_list);

    func_sum<double> ds(func_curried, x_list.size(), y_list.size() );
    std::cout<<ds.calc_section_sum(px, qx, py, qy)<<std::endl;
    return 0;
}
ABC106 D:AtCoder Express 2

問題はリンクを参照。この問題はf(x_i, y_j)  x_i = L_i, y_i = R_iのときそのペア(L_i, R_i)の存在する個数を返す関数とすれば2次元累積和で解ける。

#include <iostream>
int main()
{
    int N, M, Q;
    std::cin >> N >> M >> Q;
    std::vector<std::vector<int>> table(N, std::vector<int>(N, 0));
    for (int i = 0; i < M; i++) {
        int L, R;
        std::cin >> L >> R;
        table[L - 1][R - 1]++;
    }

    auto func = [](int i, int j,
        const std::vector<std::vector<int>>& table) { return table[i][j]; };

    std::function<int(int, int)>  func_curried = std::bind(func,
        std::placeholders::_1,
        std::placeholders::_2,
        table);

    func_sum<int> ds(func_curried, N, N);

    for (int i = 0; i < Q; i++) {
        int p, q;
        std::cin >> p >> q;
        std::cout << ds.calc_section_sum(p, q, p, q) << std::endl;
    }
    return 0;
}

しゃくとり法

用途

  • 解きたい問題が以下のいずれかの性質を満たす時、条件C(l, r)を満たすl,rの数え上げ、もしくは評価関数f(l,r)の最大化・最小化
    • l_p \leq l\leq r \leq r_qとして、C(l, r)が条件を満たすならばC(l_p, r_q)も条件を満たす(以降、これを上位区間条件と呼ぶ)
      • 例えば、整数の数列x_1, x_2, \cdots, x_nの中からその和がS以上となる連続する部分数列を数え上げる問題
    • l \leq l_p \leq r_q \leq rとして、C(l, r)が条件を満たすならばC(l_p, r_q)も条件を満たす(以降、これを部分区間条件と呼ぶ)
      • 例えば、整数の数列x_1, x_2, \cdots, x_nの中からその和がS以下となる連続する部分数列を数え上げる問題

アルゴリズム

上位区間条件を満たす問題と部分区間条件を満たす問題とで若干異なる(うまいことやれば、共通にできるかもしれない)

  1. 条件を満たさない間rをインクリメント(右端まで来てインクリメントできない場合はそのまま)
  2. 条件を満たしたところで、インクリメントを止め、その時点のl, rを記録(この時点でlp = l, r \leq r_pを満たす(l_p, r_p)は全て条件を満たす)
  3. 条件が満たさなくなるまでlをインクリメント
  4. この1~3の操作をrがインクリメントできなくなり、かつ、条件C(l, r)を満たさなくなるまで繰り返す
  1. 条件を満たす間rをインクリメント(右端まで来てインクリメントできない場合はそのまま)
  2. 条件を満たさなくなったところで、インクリメントを止め、その時点のl, rを記録(この時点でlp = l, r_p \leq rを満たす(l_p, r_p)は全て条件を満たす)
  3. 条件を満たすまでlをインクリメント
  4. この1~3の操作をrがインクリメントできなくなり、かつ、lrと同じになるまで繰り返す

抽象クラスの実装

こちらもテンプレートメソッドパターンで実装する。上述の通り上位区間条件と部分区間条件を分けて考えるので、それぞれ抽象クラスを作る。

#include <functional>
 
//ある区間が条件を満たすならば、
//それを包含する上位区間も条件を満たす問題を解くクラス
template <typename T>
class abstract_inchworm_superset
{
private:
    virtual void move_left_pointer() = 0;//lを右に動かす処理
    virtual void move_right_pointer() = 0;//rを右に動かす処理
    virtual bool is_cond() = 0; //条件成立か否かを返す
    virtual bool is_right_end() = 0;//rが右端に到達したかを返す
    virtual std::pair<T, T> no_next() = 0;//次の要素がない場合に終了値を返す
    virtual std::pair<T, T> get_LR() = 0;//lとrを返す
public:
    virtual std::pair<T, T> next(); //条件を満たすデータを返し、ポインタを次に進める
};
 
template <typename T>
std::pair<T, T> abstract_inchworm_superset<T>::next(){
    while (is_cond() == false && is_right_end() == false) {
            move_right_pointer();
    }
    if (is_cond() == false)
        return no_next();
    std::pair<T, T> ret = get_LR();
    move_left_pointer();
    return ret;
}
 
//ある区間が条件を満たすならば、
//それが包含する部分区間も条件を満たす問題を解くクラス
template <typename T>
class abstract_inchworm_subset
{
private:
    virtual void move_left_pointer() = 0;//lを右にインクリメントする処理
    virtual void move_right_pointer() = 0;//rを右にインクリメントする処理
    virtual bool is_cond() = 0; //条件成立か否かを返す
    virtual bool is_right_end() = 0;//rが右端に到達したかを返す
    virtual std::pair<T, T> no_next() = 0;//次の要素がない場合に終了値を返す
    virtual std::pair<T, T> get_LR() = 0;//lとrを返す
    virtual bool is_left_eq_right() = 0;//lとrが一致しているかを返す
public:
    virtual std::pair<T, T> next(); //条件を満たすデータを返し、ポインタを次に進める
};
 
template <typename T>
std::pair<T, T> abstract_inchworm_subset<T>::next() {
    std::pair<T, T> ret;
    while (is_right_end() == false && is_cond() == true) {
        move_right_pointer();
    }
    ret = get_LR();
    if (is_left_eq_right() == true) {
        return no_next();
    }
    
    move_left_pointer();
    return ret;
}

具体的な問題を解く実装

POJ 3061:Subsequence

蟻本に例として載っている和がS以上になる部分数列を数え上げる問題。上位区間条件型のしゃくとり法で解ける。

#include <iostream>
#include <algorithm>
#include <vector>

class subsequence :
    public abstract_inchworm_superset<int>
{
private: 
    int Left, Right;
    long long S;
    int N;
    long long LRsum;
    std::vector<long long> sequence;
    virtual void move_left_pointer();
    virtual void move_right_pointer();
    virtual bool is_cond();
    virtual bool is_right_end();
    virtual std::pair<int, int> no_next();
    virtual std::pair<int, int> get_LR();
public:

    void set_S(long long S);
    subsequence(std::vector<long long>& sequence, int N, int start);
};

subsequence::subsequence(std::vector<long long>& sequence, int N, int start):
    sequence(sequence),
    N(N),
    Left(start),
    Right(start),
    S(0),
    LRsum(0)
{
}

void subsequence::set_S(long long S) {
    this->S = S;
}

void subsequence::move_left_pointer() {
    LRsum -= sequence[Left];
    Left++;
}

void subsequence::move_right_pointer() {
    if (Right < N) {
        LRsum += sequence[Right];
        Right++;
    }
}

bool subsequence::is_cond() {
    if (LRsum >= S)
        return true;
    return false;
}

bool subsequence::is_right_end() {
    if (Right >= N)
        return true;
    return false;
}

std::pair<int, int> subsequence::no_next() {
    return std::make_pair(-1, -1);
}

std::pair<int, int> subsequence::get_LR() {
    return std::make_pair(Left, Right);
}

int main()
{

    int N;
    std::cin >> N;
    for (int k = 0; k < N; k++) {
        std::vector<long long> sequence;
        int n;
        long long S;
        std::cin >> n >> S;
        for (int i = 0; i < n; i++) {
            long long a;
            std::cin >> a;
            sequence.push_back(a);
        }

        subsequence seq(sequence, sequence.size(), 0);
        seq.set_S(S);
        int min_len = INT_MAX;
        while (1) {
            std::pair<int, int> LR = seq.next();
            if (LR.first == -1) break;
            min_len = std::min(min_len, LR.second - LR.first);
        }
        if (min_len == INT_MAX)
            std::cout << 0 << std::endl;
        else
            std::cout << min_len << std::endl;
    }
    return 0;
}
ARC098 D:Xor Sum 2

部分区間条件型のしゃくとり法で解けると気付くのが難しい問題。気づいてしまえば、普通にしゃくとり法を適用するだけ(部分区間を数え上げる問題なので、何も考えずしゃくとり法を適用してACしちゃったって人も結構いるかも)。なぜ、しゃくとり法で解けるかは解説参照。
ここで注意するのは、関数is_cond()のtrueを返す条件が現在指しているrの更に一つ先のrにおいて、条件が成立している場合としている。つまり、right pointer(r)を動かしても条件成立するなら、実際にright pointerを動かすという処理になっている。なお、この実装では区間[l, r)という半区間で扱っているため、一つ先のrrが指示しているインデックスの値を足すこと(XORを取ること)に当たる。

#include <vector>
#include <stdio.h>
class xor_sum :
    public abstract_inchworm_subset<int>
{
private:
    int Left, Right;
    
    int N;
    long long LR_XOR;
    long long LRsum;
    std::vector<long long> sequence;
 
    virtual void move_left_pointer();
    virtual void move_right_pointer();
    virtual bool is_cond();
    virtual bool is_right_end();
    virtual std::pair<int, int> no_next();
    virtual std::pair<int, int> get_LR();
    virtual bool is_left_eq_right();
public:
 
    xor_sum(std::vector<long long>& sequence, int N, int start);
};
 
xor_sum::xor_sum(std::vector<long long>& sequence, int N, int start) :
    sequence(sequence),
    N(N),
    Left(start),
    Right(start),
    LR_XOR(0),
    LRsum(0)
{
}
 
 
void xor_sum::move_left_pointer() {
    LRsum -= sequence[Left];
    LR_XOR ^= sequence[Left];
    Left++;
 
}
 
void xor_sum::move_right_pointer() {
    if (Right < N) {
        LRsum += sequence[Right];
        LR_XOR ^= sequence[Right];
        Right++;
    }
}
 
bool xor_sum::is_cond() {
    if (Left >= N) 
        return false;
    //Rightが1つ先に行っても条件を満たすか
    if ((LRsum + sequence[Right]) == (LR_XOR ^ sequence[Right]))
        return true;
    return false;
}
 
bool xor_sum::is_right_end() {
    if (Right >= N)
        return true;
    return false;
}
 
std::pair<int, int> xor_sum::no_next() {
    return std::make_pair(-1, -1);
}
 
std::pair<int, int> xor_sum::get_LR() {
    return std::make_pair(Left, Right);
}
 
bool xor_sum::is_left_eq_right() {
    return (Left == Right);
}
 
 
int main()
{
    int N;
    std::vector<long long> A;
    scanf("%d", &N);
 
    for (int i = 0; i < N; i++) {
        long long a;
        scanf("%lld", &a);
        A.push_back(a);
    }
    xor_sum xs(A, N, 0);
    long long count = 0;
    while (1) {
        std::pair<int, int> LR = xs.next();
        if (LR.first == -1) break;
        count += LR.second - LR.first ;
    }
    printf("%lld\n", count);
    return 0;
}

おわりに

今回は勉強のため、本質的なアルゴリズムを記述した抽象クラスと実際の問題を解く具象クラスを分けましたが、実際のコンテスト時はそんなことせずに、普通に書いた方が速いと思います。