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

2020/2/2追記

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

はじめに

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

yamagensakam.hatenablog.com

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

やったこと

といっても、本質的には観測モデル $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