ガウス過程の導関数を求める

はじめに

ガウス過程回帰による予測値の導関数をGPyのpredictive_gradientsメソッドによって簡単に得られることがわかったので記しておこうと思います．

ガウス過程の導関数

まずガウス過程によって推定された関数の導関数が求められることを確認します．こちらのページを参考にしました．

まず，ガウス過程回帰における予測値(事後期待値)は

で計算できます．ここで，

は新規観測値の目的変数，は新規観測値の説明変数， はデータ， はデータの説明変数と新規観測値の目的変数のカーネル値(データは複数あるのでベクトル)， はデータの説明変数のカーネル行列，はデータの目的変数を示します．そしてこれらをで微分すると，

これにより，ガウス過程回帰の事後期待値の導関数はカーネル値のベクトルの導関数を用いて求められることが分かります．

GPyによる導関数値の取得

まずはデータを作って可視化します．サインカーブにでもしてみましょう．

import GPy
import numpy as np
import matplotlib.pyplot as plt

N=100
np.random.seed(0)
x = np.random.normal(0, 1, size=N)
y = np.sin(x) + np.random.normal(0, 0.1, size=N)

kernel = GPy.kern.RBF(1)

model = GPy.models.GPRegression(x.reshape(-1, 1), y.reshape(-1, 1), kernel=kernel)
model.optimize()
model.plot()
plt.show()

GPyではガウス過程回帰の導関数がpredictive_gradientsメソッドで簡単に求められます．

このデータの導関数値を数値解析的に求めます．

x_new = np.array([[0.0]])
epsilon = 1.0e-5
x0_minus_0 = model.predict_quantiles(x_new - epsilon, quantiles=(5, 50, 95))[1][0][0]
x0_plus_0 = model.predict_quantiles(x_new + epsilon, quantiles=(5, 50, 95))[1][0][0]
print((x0_plus_0 - x0_minus_0) / (2 * epsilon))
# 0.9633798891961475

ここでは小さい値だけ前後にずらした2点の傾きを求めています． 今度はpredictive_gradientsメソッドを使ってみます．

print(model.predictive_gradients(x_new)[0][0][0][0]) 
# 0.9633798901678498

ちゃんと求められていますね(数値解析的に求めたものを答えとするのってどうなんだろう…)．GPyのgithubからも分かるように予測値の分散の導関数(導関数の分散ではない)も一つ目の添字を1にすることで得られます．

まとめ

今回はガウス過程回帰の予測値の導関数を求められることの確認とGPyライブラリによる計算方法の紹介でした．