前回 の続き。GaussianProcesses.jlを使ってガウス過程回帰のハイパーパラメータ推定を行ってみる。 「ガウス過程と機械学習」の図3.20(a)を例とする。

xxs = [-0.5, 0.5, 1, 1.4, 3]
yys = [0.7, 1.8, 1.7, 2.3, 1]

scatter(xxs, yys)
0 1 2 3 0.9 1.2 1.5 1.8 2.1 y1

平均0のガウスカーネルのガウス過程に(ノイズ項込みで)フィッティングさせるのは非常に簡単で、optimize! を呼べば良いだけ。

gp = GP(xxs, yys, MeanZero(), SE(0.0, 0.0))
optimize!(gp)

Results of Optimization Algorithm
 * Algorithm: L-BFGS
 * Starting Point: [-2.0,0.0,0.0]
 * Minimizer: [-1.250611250087999,0.5939445695924219, ...]
 * Minimum: 5.464078e+00
 * Iterations: 10
 * Convergence: true
   * |x - x'| ≤ 0.0e+00: false 
     |x - x'| = 1.80e-09 
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: false
     |f(x) - f(x')| = 1.63e-16 |f(x)|
   * |g(x)| ≤ 1.0e-08: true 
     |g(x)| = 1.23e-14 
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 26
 * Gradient Calls: 26

パラメーターを確認すると、

print(gp)

GP Exact object:
 Dim = 1
 Number of observations = 5
 Mean function:
   Type: MeanZero, Params: Float64[]
 Kernel:
   Type: SEIso{Float64}, Params: [0.593945, 0.233886]
 Input observations = 
[-0.5 0.5 1.0 1.4 3.0]
 Output observations = [0.7, 1.8, 1.7, 2.3, 1.0]
 Variance of observation noise = 0.08198471101193598
 Marginal Log-Likelihood = -5.464

「ガウス過程と機械学習」では、カーネルは $$ k(\mathbf{x}, \mathbf{x}^\prime | \bm{\theta}) = \theta_1 \exp \left( - \frac{|\mathbf{x} - \mathbf{x}^\prime|^2}{\theta_2}\right) + \theta_3 \delta(\mathbf{x}, \mathbf{x}^\prime) $$ となっていて、推定結果は \((\theta_1, \theta_2, \theta_3)=(1.596,6.560,0.082)\) とされている。 \(\theta_3\)はVariance of observation noiseの0.08198471101193598とそのまま一致。 SEIsoのカーネル関数は $$ k(\mathbf{x}, \mathbf{x}^\prime) = \sigma^2 \exp \left( - \frac{(\mathbf{x} - \mathbf{x}^\prime)^\top(\mathbf{x} - \mathbf{x}^\prime))}{2l^2}\right) $$ で、パラメーターは \(l\) と \(\sigma\) の対数をとった ll と lσ だから \(\theta_1 = \exp(2 * l\sigma), \theta_2=2 * \exp(2 * ll)\) という変換で対応できる。 SEIso のメンバー変数 σ2 と ℓ2 を使って確認すると、

println("theta1: ", gp.kernel.σ2)
println("theta2: ", gp.kernel.ℓ2 * 2)

theta1: 1.5964347486496255
theta2: 6.560299908989858

となり、同じ結果が得られている。

plot(gp; xlabel="gp.x", ylabel="gp.y", legend=false, fmt=:png)
-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 gp.x gp.y

上だと点が端ギリギリに配置され結構見ずらい(+ガウス誤差が含まれていない？)ので、ガウス観測誤差も含めた予測結果を自分でプロットすると、

xtest = reshape(collect(range(-1, stop = 3.5, length = 100)), 1, :)
μ, σ2 = predict_y(gp, xtest)
dists = Normal.(μ, sqrt.(σ2))

p = [0.025, 0.975]
qt = hcat(map(x -> quantile.(x, p), dists)...)

plot(xtest[:], qt[1, :], fillrange = qt[2, :], fillalpha = 0.3,
    label = "", linewidth = 0)
scatter!(xxs, yys, label = "")
plot!(xtest[:], mean.(dists), linewidth = 2, 
      linestyle = :dash, label = "")
-1 0 1 2 3 -0.5 0.0 0.5 1.0 1.5 2.0 2.5

ノートブック: https://github.com/matsueushi/notebook_blog/blob/master/gaussianprocess-julia.ipynb