ガウス過程回帰を実装して遊んだ

ガウス過程の本が発売されるということで復習として簡単にメモしようとおもったけど、はてなのTeX記法にうまく馴染めなかった¹ので数式は載せずに結果だけ簡単に貼り付けて日記としておきます。

結果

ガウス過程回帰、ハイパーパラメタをうまく決めないとゴミみたな結果が得られるようです。適当な回数イテレーションを回してハイパーパラメタを決めましょう。やり方はいろいろあると思いますが、MCMCで決めるのがいいかなと考えています。勾配法はうまくいったためしがないので。。。

f:id:kenbell1988:20190211143121p:plain — ガウス過程回帰の結果.(左)入力、(中)適当なパラメタ、(右)よさそうなパラメタ

f:id:kenbell1988:20190211143126p:plain — ガウス過程回帰の結果.(左)入力、(中)適当なパラメタ、(右)よさそうなパラメタ

f:id:kenbell1988:20190211143816p:plain — 2次元のシグモイド関数でやってみた結果。右上と左下はうまくいってません。右下は良さそうに見えます。

実装

Julia langで実装しました。でももしかしたらあんまりJuliaっぽくない書き方かもしれません。
３次元くらいまでは可視化して人間が理解できるのでやってみました。ここでは2次元の結果までのせておきます。

どっかで自分で使うこともあるとおもうので、コードはgithubに放り投げておきました。
github.com

下のコードは1次元のやつ。

using PyCall, PyPlot
using Distributions

function kronecker(i, j)
    if i == j
        return 1
    else
        return 0
    end
end

function kernel(x1, x2, theta)
    # 任意の字数のベクトルの写像
    sigf = theta[1]
    sigls = theta[2:end]
    L    = inv(Diagonal(sigls))^2
    coef = ((x1 - x2)' * L * (x1 - x2))[1] # 1x1になってしまうため、要素をとりだして数値にする
    return (sigf^2) * exp(-1/2.0 * coef)
end

function gram_matrix(x1, x2, theta)
    n = size(x1)[1]
    m = size(x2)[1]
    K = zeros(n, m)    for i=1:n
        for j=1:m
            K[i, j] = kernel(x1[i], x2[j], theta) 
        end
    end
    return K
end

function log_likelihood(x, y, b, theta)
    t = y
    n = size(x)[1]
    K= gram_matrix(x, x, theta)
    C = K + eye(n)/b
    tCt = t' * inv(C) * t
    logdetC = logdet(C)  # log(det(C))
    return  -1/2.0 * tCt - 1/2.0 * logdetC - n/2.0 * log(2*pi)
end

function regression(x_obs, y_obs, b, th)

    # 共分散行列の計算
    K = gram_matrix(x_obs, x_obs, th)
    n = size(y_obs)[1]
    C_N = K + eye(n)/b

    # 平均値の計算
    m_list = []
    for j = 1:length(x)
        k = [kernel(x_obs[i, :], x[j], th) for i =1:length(y_obs)]
        m =   k' *  inv(C_N)  * y_obs
        push!(m_list, m[1])
    end

    # 分散の計算
    sig2_list = []
    for j = 1:length(x)
        k = [kernel(x_obs[i, :], x[j], th) for i =1:length(y_obs)]
        c = kernel(x[j], x[j], th) + 1/b
        push!(sig2_list, c - k' * inv(C_N) * k)
    end
    return m_list, sig2_list 
end

function y(x)
    return cos(pi/10*x) + exp(x/20) - 0.5 * (x/20)^4
end

#　適当にデータをつくる
n = 40
b = 20 
x_obs = rand(n)*50
sig = 1
y_obs = y.(x_obs) + rand(Normal(0, sig), n);
x = linspace(0, 50, 101);

# 入力の可視化
figure("result")
clf()
plot(x, y.(x), c="springgreen", label="ground truth")
plot(x_obs, y_obs, "o", c="black", label="observe")
xlabel("x")
ylabel("y")
title("Input data")
legend()
show()

# 適当な値
b = 1.0          # うまくいかない
th = [1.0, 1.0]　# うまくいかない
m_list, sig2_list = regression(x_obs, y_obs, b, th)

figure("result_1")
clf()
plot(x, y.(x), c="springgreen", label="ground truth")
plot(x, m_list, "--", label="mean")
plot(x, m_list - 3*sig2_list, "--", c="orange")
plot(x, m_list + 3*sig2_list, "--", c="orange", label="3σ")
plot(x_obs, y_obs, "o", c="k", label="observe")
title("Random params")
legend()
show()

# 最適化した値
b = 1.1411959945407517
th = [2.388171072439165, 6.052294424769037]
m_list, sig2_list = regression(x_obs, y_obs, b, th)

figure("result_2")
clf()
plot(x, y.(x), c="springgreen", label="ground truth")
plot(x, m_list, "--", label="mean")
plot(x, m_list - 3*sig2_list, "--", c="orange")
plot(x, m_list + 3*sig2_list, "--", c="orange", label="3σ")
plot(x_obs, y_obs, "o", c="k", label="observe")
title("Optimized parameters")
legend()
show()