初っ端ですがGauss過程回帰を実装してみます。今回は全面的にこちらの書籍を参考にしていきます。

-『ガウス過程と機械学習』持橋大地・大羽成征 MLP

作者:持橋大地,大羽成征
発売日: 2019/03/09
メディア: 単行本（ソフトカバー）

1. 理論
- 1.1. Gauss過程とは
- 1.2. Gauss過程回帰
2. 実装
- 2.1. Markov Chain Monte Carlo (MCMC)を用いたハイパーパラメータの決定
- 2.2. GPR実装
3. 結果
- 3.1. 回帰結果
- 3.2. MCMCの様子
4. おまけ (MATLAB ソースコード)
5. 参考文献
6. さいごに・感想

1. 理論

理論に関しては初学者の自分がへたに解説する間でもないと思いますので、ざっくりといきます。何か間違っていたらコメント欄で優しく教えてください。

1.1. Gauss過程とは

いま観測 $\textbf{y} = (y_1, y_2, \cdots, y_N) ^ T$ を $\textbf{x}_1, \textbf{x}_2, \cdots, \textbf{x}_N$ を用いて回帰することを考えます。

$\textbf{y}$ の予測を $\hat{\textbf{y}} =(\hat{y}_1, \hat{y}_2, \cdots, \hat{y}_N)$ とおき、 $\textbf{x}$ の基底を並べた特徴ベクトル を $\boldsymbol{\phi} ( \textbf{x} ) = (\phi_0(\textbf{x}), \phi_1(\textbf{x}), \phi_2(\textbf{x}), \cdots, \phi _H(\textbf{x}))^T$ , 重みベクトルを $\textbf{w} =(w_1, w_2, \cdots, w_H)$ とおくとこのようにかけます。

$\displaystyle \begin{pmatrix} \hat{y}_{1} \\ \hat{y}_{2} \\ \vdots \\ \hat{y}_{N} \end{pmatrix} = \begin{pmatrix} 1 & \phi_1 (\textbf{x}_1) & \phi_2 (\textbf{x}_1) & \cdots & \phi_H (\textbf{x}_1) \\ 1 & \phi_1 (\textbf{x}_2) & \phi_2 (\textbf{x}_2) & \cdots & \phi_H (\textbf{x}_2) \\ \vdots & & & & \vdots \\ 1 & \phi_1 (\textbf{x}_N) & \phi_2 (\textbf{x}_N) & \cdots & \phi_H (\textbf{x}_N) \\ \end{pmatrix} \begin{pmatrix} w_{1} \\ w_{2} \\ \vdots \\ \vdots \\ w_{H} \end{pmatrix}$

ベクトルで書くと

$\displaystyle \hat{\textbf{y}} = \boldsymbol{\Phi} \textbf{w}$

例えば $\phi (\textbf{x})$ として多項式や動径基底関数などを考えるわけですが、例えば動径基底函数に対応する $\textbf{w}$ は次元の呪いによって $H^D$ のオーダーで指数的に増えていってしまい、簡単に計算量が爆発します( $\boldsymbol{\Phi}$ の基底関数 $\phi$ が横に指数的にどんどん長くなっていく)。

これはあまり嬉しくないので、重みが $\displaystyle \textbf{w} \sim \mathcal{N}(\textbf{0}, \lambda ^2 \textbf{I})$ のような確率分布から生成されたとして、これ積分消去できたらいいなと考えます。いま単純のため $\hat{\textbf{y}} = \textbf{y}$ とすると $\textbf{y}$ の共分散 $\Sigma$ は

$\begin{eqnarray} \Sigma &=& E[\textbf{y} \textbf{y}^T] - E[\textbf{y}] E[\textbf{y}] ^T \\ &=& E[(\boldsymbol{\Phi}\textbf{w})(\boldsymbol{\Phi}\textbf{w}) ^T] \\ &=& \boldsymbol{\Phi} E[\textbf{w}\textbf{w} ^T] \boldsymbol{\Phi} ^T = \lambda ^2 \boldsymbol{\Phi} \boldsymbol{\Phi} ^T \end{eqnarray}$

つまり $\textbf{y} \sim \mathcal{N} (\textbf{0}, \lambda ^2 \boldsymbol{\Phi} \boldsymbol{\Phi} ^T)$ となって $\textbf{y}$ は平均0の多変量Gauss分布に従います。そして狙った通り $\textbf{w}$ は消えてくれ、高次元の計算をわざわざする必要はなくなりました。

この確率過程 $\textbf{y}$ をGauss過程 (Gaussian process)と呼びます。正確な定義は以下のようになります。

Gauss過程
あらゆる自然数 $N$ で、入力 $\textbf{x}\_1, \textbf{x}\_2, \cdots, \textbf{x}\_N$ に対応する出力

$\displaystyle \textbf{f} = ( f(\textbf{x}_1), f(\textbf{x}_2), \cdots, f(\textbf{x}_N))$

が平均 $\displaystyle \boldsymbol{\mu} = (\mu(\textbf{x}_1), \mu(\textbf{x}_2), \cdots, \mu(\textbf{x}_N)$ ), 共分散行列 $\textbf{K}$ とする多変数Gauss分布に従うとき、これをGauss過程(Gaussian process) と呼ぶ。
ただし $\textbf{K}$ は $i, j$ 要素に関して以下が成り立つ行列である( カーネル行列(kernel matrix)あるいはグラム行列(Gram matrix)という)。

$\displaystyle (K)_{ij} = k(\textbf{x}_i, \textbf{x}_j) = \boldsymbol{\phi} (\textbf{x}_i) ^T \boldsymbol{\phi} (\textbf{x}_j)$

このGram行列というやつが重要です。上の例では $\textbf{y} \sim \mathcal{N} (\textbf{0}, \lambda ^2 \boldsymbol{\Phi} \boldsymbol{\Phi} ^T)$ の定数倍を無視して $\boldsymbol{\Phi} \boldsymbol{\Phi} ^T$ の部分を $\textbf{K}$ として議論します。

Gram matrixがすごいのは、kernelの式さえうまく指定すれば特徴ベクトル(基底ベクトル)の計算をいちいちしないでよくなるという点です(詳細はMLP ガウス本を参照してください)。これをカーネルトリック (Kernel trick)といい、須くこのtrickを用いる手法をkernel法といいます。

Kernelの式は勿論なんでもいいわけではなく、「うまく」指定しなくてはいけませんが、これは $\textbf{K}$ がGauss分布の共分散として機能するようにしてくださいねということで、つまり $\textbf{K}$ を計算したとき対称かつ正定値にならなくてはなりません。

Gauss過程に用いるkernelにはいろいろと種類がありますが、一番オーソドックスなRBF kernel (radial basis function kernel, Gaussian kernel)を示します。

$\displaystyle k(\textbf{x}, \textbf{x'}) = \theta_1 \exp \left( - \frac{\left| \textbf{x} - \textbf{x'} \right| ^2} {\theta_2} \right)$

RBF kernelは、 $\textbf{x}, \textbf{x'}$ の値が近かったら必然的に値が大きくなります。つまり「入力が近かったら出力も近くなる」という「類似度」の意味が込められているというわけですね。このほかにも指数カーネル(exponential kernel)、周期カーネル(periodic kernel)などがあります。

Kernelに関しての議論は難しいらしいのでこれ以上あまり深入りしないことにしますが(再生核Hilbert空間での内積とか)、これらのkernelより作られるGram matrixは正定値であることが保証されています(証明はあんまり理解していない)。よく考えたらすごいですね。そしてKernelどうしのの和・積なども半正定値性を持ちます。これも普通にすごいことなので証明をきちんと理解したいです......

さて、もう一度 $\textbf{y} \sim \mathcal{N} (\textbf{0}, \lambda ^2 K)$ という式に立ち返ってみます。この式のお気持ちとしては、入力 $\textbf{x}_1, \textbf{x}_2, \cdots, \textbf{x}_N$ を決めると、出力 $\textbf{y}$ の各次元の予測は $-\infty$ から $\infty$ まで全ての実数値の個数次元(=無限次元)のGauss分布からサンプリングされるということです。函数の値が無限次元のGauss分布から出てくるなんてちょっと面白い考え方ですね。

1.2. Gauss過程回帰

Gauss過程を回帰問題に適応するにはどうしたらよいでしょう。ということですが、上で見た入力 $\textbf{x}_1, \textbf{x}_2, \cdots, \textbf{x}_N$ とそれぞれの観測 $\textbf{y} = (y_1, y_2, \cdots, y_N) ^ T$ を $\textbf{x}_1, \textbf{x}_2, \cdots, \textbf{x}_N$ があるとき、新たな入力 $\textbf{x}^{*}$ とその予測値 $y^{*}$ を考えます。

すると、もとの観測と新しい予測を合わせた $\textbf{y'} = (y_1, y_2, \cdots, y_N, y^{*}) ^ T$ は、もとの入力と新しい入力を合わせた $\textbf{x}_1, \textbf{x}_2, \cdots, \textbf{x}_N, \textbf{x}^{*}$ から計算できるGram matrix $\textbf{K'}$ ( $(N+1)$ 次正方行列)を共分散行列に持つGauss過程です。新しい入力を付け加えてもGauss過程というのがミソですね。

$\displaystyle \begin{pmatrix} \textbf{y}_1 \\ \textbf{y}_2 \end{pmatrix} \sim \left( \begin{pmatrix} \boldsymbol{\mu}_1 \\ \boldsymbol{\mu}_2 \end{pmatrix} , \begin{pmatrix} \boldsymbol{\Sigma}_{11} & \boldsymbol{\Sigma}_{12} \\ \boldsymbol{\Sigma}_{21} & \boldsymbol{\Sigma}_{22} \end{pmatrix} \right)$

このとき

$\displaystyle p(\textbf{y}_2|\textbf{y}_1) = \mathcal{N} \left( \boldsymbol{\Sigma}_{21} \boldsymbol{\Sigma}_{11}^{-1} \textbf{y}_1, \boldsymbol{\Sigma}_{22} - \boldsymbol{\Sigma}_{21} \boldsymbol{\Sigma}_{11}^{-1} \boldsymbol{\Sigma}_{12} \right)$

が成り立つので(証明はMLP本に加えPRMLのはじめの方にも載っている)、Gauss過程での予測は

$\begin{eqnarray} \textbf{k}_{*} &=& \left( k(\textbf{x}^{*}, \textbf{x}_1), k(\textbf{x}^{*}, \textbf{x}_2), \cdots, k(\textbf{x}^{*}, \textbf{x}_N) \right)^T \\ k_{**} &=& k(\textbf{x}^{*}, \textbf{x}^{*}) \end{eqnarray}$

として

$p(y^{*}|\textbf{x}^{*}, \mathcal{D}) = \mathcal{N} \left( \textbf{k}_{*}^T \textbf{K}^{-1} \textbf{y}, k_{**} - \textbf{k}_{*}^T \textbf{K}^{-1} \textbf{k}_{*} \right)$

同様に、 $y$ がベクトルの場合も

$p(\textbf{y}^{*}|\textbf{x}^{*}, \mathcal{D}) = \mathcal{N} \left( \textbf{k}_{*}^T \textbf{K}^{-1} \textbf{y}, k_{**} - \textbf{k}_{*}^T \textbf{K}^{-1} \textbf{k}_{*} \right)$

要は、「インプットと観測データが与えられたとき、インプットと観測の間の函数を、どの程度の信頼がおけるか、という情報も込めて予測できる」ということなので、どの程度の精度をもって予測できるかわかるのは非常に嬉しいですよね。

Kernel函数ですが、今回の実験にはこちらのkernelを使用します( $\delta$ はKronecker's deltaで、Gram matrixの対角項に加える誤差項を意味します)。

$\displaystyle k(\textbf{x}, \textbf{x'}) = \theta_1 \exp \left( - \frac{\left| \textbf{x} - \textbf{x'} \right| ^2} {\theta_2} \right) + \theta_3 \delta (\textbf{x}, \textbf{x'})$

意味としては、いま簡単のため $\displaystyle \textbf{y} = \hat{\textbf{y}}$ と仮定してやっていましたが、きちんと誤差項を入れた $\textbf{y} = \textbf{f} (\textbf{x}) + \boldsymbol{\epsilon}$ を考えると、共分散行列の対角項に誤差の分散が足され、上が導けます(詳しくはMLP本を参照)。

なのでKronecker's deltaの部分 $\delta (\textbf{x}, \textbf{x'})$ は $\textbf{x}, \textbf{x'}$ の値が等しければ1というわけではなく対角項が1と返ってくるだけです。

2. 実装

Juliaを使って実装します。

(これからこのブログではJuliaをメインで使って実装していこうと思います。メインストリームと言えばPythonでしょうが、個人的に少し書きにくいのと、自分がひねくれている性格ので一番人気のあるやつはやめときたいというだけです。)

データセットとしてはMLP本のサポートページからダウンロードできるデータセットに加えて、自分で適当にデータを作ってやってみます。

$\displaystyle y = A \sin \left( \frac{2\pi (x-x_0)}{T} \right) + B \epsilon$

$(A, B, x_0 , T) = (50, 25/2, -15, 120)$ として、 $\epsilon \sim \mathcal{N}(0,1)$ としてエラーを加えて得た200サンプルです。プロットはこんな感じです。

f:id:gvsk:20200306231401p:plain — fig. 2-1

そして、せっかくならなにか実データでやってみたいということで、少し前に「ベテルギウスBetelgeuseが暗くなっている」というニュースを読んで気になっていたので、データをAmerican Association of Variable Star Observers(AAVSO)からダウンロードしてやってみました。

2.1. Markov Chain Monte Carlo (MCMC)を用いたハイパーパラメータの決定

いざ実装！と行きたいわけですが、その前にやっておかなくてはいけないことがあります。これがハイパーパラメータ(hyperparameters)の推定です。ガウス過程での「学習」はkernelの式の中のハイパーパラメータを通してのみ行われます。

$\displaystyle k(\textbf{x}, \textbf{x'}) = \theta_1 \exp \left( - \frac{\left| \textbf{x} - \textbf{x'} \right| ^2} {\theta_2} \right) + \theta_3 \delta (\textbf{x}, \textbf{x'})$

この式の中の $\theta_1, \theta_2, \theta_3$ ですね。ハイパーパラメータの値を調整しないとどうなるでしょうか。結果を一部先取りしてしまいますが、 $(\theta_1, \theta_2, \theta_3) = (0.1, 0.002, 0.003)$ としてやってみます。

f:id:gvsk:20200306224359p:plain — fig. 3-1

方法自体は正しいにもかかわらず、over-fittingしてイケてない回帰になっていますね。というわけで、ハイパーパラメータは是非いい感じにしたい。

このためにパラメータから尤度を計算し、尤度を最大化する方法をとります。

$\displaystyle \begin{eqnarray} p(\textbf{y}|\textbf{X}, \boldsymbol{\theta}) &=& \mathcal{N} (\textbf{y}|\textbf{0}, \textbf{K}_\boldsymbol{\theta} ) \\ &=& \frac{1}{(2\pi)^{N/2}} \frac{1}{\left| \textbf{K}_\boldsymbol{\theta} \right| ^{1/2}} \exp(-\frac{1}{2} \textbf{y}^T \textbf{K}_\boldsymbol{\theta} ^{-1} \textbf{y}) \end{eqnarray}$

よって $\log$ をとって対数尤度にすれば

$\displaystyle \begin{eqnarray} \log p(\textbf{y}|\textbf{X}, \boldsymbol{\theta}) &=& - \frac{N}{2} \log (2\pi) - \frac{1}{2} \log \left| \textbf{K}_\boldsymbol{\theta} \right| - \frac{1}{2} \textbf{y}^T \textbf{K}_\boldsymbol{\theta} ^{-1} \textbf{y} \\ &\propto& - \log \left| \textbf{K}_\boldsymbol{\theta} \right| - \textbf{y}^T \textbf{K}_\boldsymbol{\theta} ^{-1} \textbf{y} + {\rm const.} \end{eqnarray}$

上記を最大化すればいいわけです。よく用いられる手法としてMarkov chain Monte Carlo (MCMC)や、勾配法 gradient methodがあります。

MLP本では、この式の勾配が解析的に求められるので、後者の解説がされていました。しかしながら今回は、勾配法は局所解に陥りやすいが、MCMCならそれを避けながら、解析的に求められない問題にも適用可能だということでこちらで実装してみます。本当は、MCMCを以前授業でやったことがあっただけです(MCMCだけでも記事が書けるのでいずれ書こうと思います)

MCMCの詳しい説明は今回踏み込みませんが、サンプリングの難しい確率分布からMarkov連鎖を用いてうまいことサンプリングする手法です。今回は最もクラシカルなMetropolis-Hastings法を用いて試みます。本当はMarkov連鎖が収束するとか詳細釣り合い条件とかいろいろ議論しないといけないのですが、手順だけ書いておきます。

Metropolis Hastings法

初期値 $\textbf{x}^{(0)}$ を決定する。
について以下を繰り返す。
1. $\textbf{y}$ を $q(\textbf{y}|\textbf{x})$ より発生させる。ここでは
  $\displaystyle \textbf{y}=\textbf{x}+\boldsymbol{\epsilon}, ~~~\boldsymbol{\epsilon} \sim \mathcal{N} (\textbf{0},\,\sigma^{2} \textbf{I})$
  に従って出していく(random walk) 。
2. 発生させた $\textbf{y}$ に対して、採択確率 $\alpha (\textbf{x}, \textbf{y}) = \min\{ 1, \frac{P(\textbf{y})q(\textbf{x}|\textbf{y})}{P(\textbf{x})q(\textbf{y}|\textbf{x})} \}$ で採択し、それ以外の場合は棄却して $\textbf{x}$ を再度採用する。

2.2. GPR実装

回帰の実装自体はそんなに難しくありません。

## Load data
using LinearAlgebra
using DelimitedFiles
using Distributions
using Plots
gr()
data = readdlm("gpr.dat",'\t',Float64);
N = 5; #number pf training data points
A = data[1:N,:];
xtrain = A[:,1]; ytrain = A[:,2];

## Define Gaussian Kernel
function GaussianKernelMat(x1, x2, θ1, θ2)
    N1 = size(x1,1); N2 = size(x2,1);
    x_d = x1 * ones(1,N2) - (x2 * ones(1,N1))';
    x_log = - x_d .^2 / θ2;
    return θ1 * (exp(1) * ones(N1,N2)) .^ x_log;
end


## Detemine hyperparameters via MCMC (Metropolis-Hastings method)
function LogLikelihood(x,y,θ1,θ2,θ3)
    Nx = length(x);
    K_θ = GaussianKernelMat(x, x, θ1, θ2) + θ3 * Matrix(1.0I,Nx,Nx);
    return LLH = - logdet(K_θ) - y' * inv(K_θ) * y;
end

step_n  = 100000; #step number of MCMC
burn_in = 5000;   #steps at the beginning to discard
θ1_0    = 1.0; θ2_0 = 1.0; θ3_0 = 0.2; #initial values
ϕ_ts    = [log(θ1_0) log(θ2_0) log(θ3_0); zeros(step_n - 1,3)]; #ϕ represents logθ
L_ts    = [LogLikelihood(xtrain,ytrain,θ1_0,θ2_0,θ3_0); zeros(step_n -1 ,1)];
step_w  = 1.0;    #step width of one MH step

for t = 2:step_n
    ϕ_tmp = ϕ_ts[t-1,:] + step_w * randn(3,1);
    θ_tmp = (exp(1) * ones(Float64,3)) .^ ϕ_tmp;
    K = GaussianKernelMat(xtrain, xtrain, θ_tmp[1],θ_tmp[2]) + θ_tmp[3] * Matrix(1.0I,N,N);

    r = rand(Float64);
    ratio = exp(LogLikelihood(xtrain,ytrain,θ_tmp[1],θ_tmp[2],θ_tmp[3]) - LogLikelihood(xtrain,ytrain,exp(ϕ_ts[t-1,1]),exp(ϕ_ts[t-1,2]),exp(ϕ_ts[t-1,3])));
    if ratio > r
        ϕ_ts[t,:] = ϕ_tmp;
        L_ts[t] = LogLikelihood(xtrain,ytrain,θ_tmp[1],θ_tmp[2],θ_tmp[3]);
    else
        ϕ_ts[t,:] = ϕ_ts[t-1,:];
        L_ts[t] = L_ts[t-1];
    end
end

#Cut burn-in steps
ϕ_ts = ϕ_ts[burn_in+1:end,:];
L_ts = L_ts[burn_in+1:end,:];

#take out maximal values
ff = findmax(L_ts); ci = ff[2];
argmax_ϕ = ϕ_ts[ci[1],:];  #ϕ values at the maximal likelihood
argmax_θ = (exp(1) * ones(Float64,3)) .^ argmax_ϕ; #θ values at the maximal likelihood


#range of input
xtest = -1:0.01:3.5;
xlen = length(xtest);

#estimate mean of y
K = GaussianKernelMat(xtrain,xtrain,argmax_θ[1],argmax_θ[2]) + argmax_θ[3] * Matrix(1.0I,N,N);
k_star = GaussianKernelMat(xtrain,xtest,argmax_θ[1],argmax_θ[2]);
y_est = k_star' * inv(K) * ytrain;

k_double_star = argmax_θ[1] + argmax_θ[3];

#variance of y
y_var = k_double_star * ones(xlen,1) - diag(k_star' * inv(K) * k_star);
y_sig = y_var .^(1/2);

#lastly, figure plot
plot(xtest, y_est, fillrange=[y_est-2y_sig], fillalpha=0.3, c=:orange)
plot!(xtest, y_est, fillrange=[y_est+2y_sig], fillalpha=0.3, c=:orange)
plot!(xtrain, ytrain, seriestype=:scatter)

3. 結果

3.1. 回帰結果

はじめに、MLP ガウス本のデータを使って実験です。全データ10点のうちはじめの5点のみ使って回帰した結果は以下のようになり、MLP ガウス本p.91の図3.20(a)を再現することができました。ハイパーパラメータの値は $(\theta_1, \theta_2, \theta_3) = (1.480, 5.832, 0.086)$ で、これも同本に記載のある勾配法を用いた値( $(\theta_1, \theta_2, \theta_3) = (1.596, 6.560, 0.082)$ )と少しずれがあるものの大方一致しました(MCMCのステップ幅をせばめればもう少しこの値に近くかもしれません)。

f:id:gvsk:20200305232213p:plain — fig. 2

次に10点すべてを用いて行なった回帰を示します。MLP本p.91の図3.20(b)にあたります。ハイパーパラメータは $(\theta_1, \theta_2, \theta_3) = (1.496, 0.6913, 0.0690)$ で、これも本の値とほぼ一致しました。複雑な回帰がいい感じでできていることがわかります。

f:id:gvsk:20200305232549p:plain — fig. 3

はじめにつくった自作のサンプルデータでも同じことをやりました。

f:id:gvsk:20200306102607p:plain — fig3-4

MCMCのおかげでハイパーパラメータも推定でき、データ点が多くても精度良く予測してくれます。データ点の近くでは分散が小さく自信のある様子を、すこし外れたところでは分散が大きく自信のない程度も教えてくれています。

最後に、Betelgeusでの実験はこの通りです。2010年から直近10年の結果を使用し、1年のうち2ヶ月分の平均を6回とって使っています。

f:id:gvsk:20200306232916p:plain — fig. 3-2

滑らかな関数では回帰できており、dimmingの様子もわかりますが、いまいちなんとも言えないですねぇ。まぁ、平均0のGauss過程でGauss kernelをそのまま使っているのでしょうがないでしょう。もう少し平均やkernelに工夫をすればマシになるかもしれません。

3.2. MCMCの様子

今回は3パラメータありますので、パラメータがMH法で動いた様子を3D scatter plotで表現しました。MLP本のデータ10点で学習したときのものです。

f:id:gvsk:20200305233912p:plain — fig. 3-3

確率分布の山の高そうなところに点が集まっているのが分かります。

4. おまけ (MATLAB ソースコード)

研究室ではいつもMATLABを使っているので、自分としてはこちらの方が書きやすいです。今回Juliaはほぼ初めてでしたが、バグ取り中にエラーの場所を確認する意味でMATLABに移植して同じことをしていたため、副産物的にMATLABコードもできてしまいました(Juliaでやる意味)。こちらも載せるだけ載せておきます。Built inの函数などを使って書けばもっと少なくかけると思います。

dir = which('sample_data_200.txt');
data = readmatrix(dir);

%load('~/Desktop/machine_learning/gpbook_data/gpr.dat')
%data = gpr;
N = 100; A = data(1:N,:);
xtrain = A(:,1); ytrain = A(:,2);

t1_0 = 1.0; t2_0 = 1.0; t3_0 = 0.2;

%% Metropolis Hastings for hyperparameter determination
MCMCsteps = 100000; burn_in = 5000;
p_ts = [log(t1_0) log(t2_0) log(t3_0); zeros(MCMCsteps - 1,3)];
L_ts = [LogLikelihood(xtrain,ytrain,t1_0,t2_0,t3_0); zeros(MCMCsteps -1 ,1)];
step_w = 0.5;

for t = 2:MCMCsteps
    p_tmp = p_ts(t-1,:) + step_w * randn(1,3);
    t_tmp = (exp(1) * ones(1,3)) .^ p_tmp;
    K = GaussianKernelMat(xtrain, xtrain, t_tmp(1),t_tmp(2)) + t_tmp(3) * eye(N);
    while any(eig(K) <= 0) % limit parameters for K to be positive definite
        p_tmp = p_ts(t-1,:) + step_w * randn(1,3);
        t_tmp = (exp(1) * ones(1,3)) .^ p_tmp;
        K = GaussianKernelMat(xtrain, xtrain, t_tmp(1),t_tmp(2)) + t_tmp(3) * eye(N);
    end
    r = rand;
    ratio = exp(LogLikelihood(xtrain,ytrain,t_tmp(1),t_tmp(2),t_tmp(3))...
        - LogLikelihood(xtrain,ytrain,exp(p_ts(t-1,1)),exp(p_ts(t-1,2)),exp(p_ts(t-1,3))));
    if ratio > r
        p_ts(t,:) = p_tmp;
        L_ts(t) = LogLikelihood(xtrain,ytrain,t_tmp(1),t_tmp(2),t_tmp(3));
    else
        p_ts(t,:) = p_ts(t-1,:);
        L_ts(t) = L_ts(t-1);
    end
end

p_ts = p_ts(burn_in+1:end,:); L_ts = L_ts(burn_in+1:end,:);
[maxL, ind] = max(L_ts);
argmax_t = (exp(1) * ones(1,3)) .^ p_ts(ind,:);


%% Draw figure
xtest = [-20:0.1:200]';
xlen = length(xtest);

K = GaussianKernelMat(xtrain,xtrain,argmax_t(1),argmax_t(2)) + argmax_t(3) * eye(N);
k_star = GaussianKernelMat(xtrain,xtest,argmax_t(1),argmax_t(2));
y_est = k_star' * inv(K) * ytrain;

k_double_star = argmax_t(1) + argmax_t(3);

y_var = k_double_star * ones(xlen,1) - diag(k_star' * inv(K) * k_star);
y_sig = y_var .^(1/2);

figure;
h = plot(xtest,y_est);
hold on
scatter(xtrain,ytrain,'x');
ar = area(xtest,[y_est - 2*y_sig, 2*y_sig, 2*y_sig]);
hold off
set(h,'Color',[0.8500 0.3250 0.0980], 'LineWidth',1.0);
set(ar(1),'FaceColor',[1 1 1],'LineStyle','None');
set(ar(2),'FaceColor',[0.8500 0.3250 0.0980],'LineStyle','None','FaceAlpha',0.25);
set(ar(3),'FaceColor',[0.8500 0.3250 0.0980],'LineStyle','None','FaceAlpha',0.25);

% plot(xtest, y_est, fillrange=[y_est-2y_sig], fillalpha=0.3, c=:orange)
% plot!(xtest, y_est, fillrange=[y_est+2y_sig], fillalpha=0.3, c=:orange)
% plot!(xtrain, ytrain, seriestype=:scatter)

%% Define functions
function GK = GaussianKernelMat(x1, x2, t1, t2)
    N1 = size(x1,1); N2 = size(x2,1);
    x_d = x1 * ones(1,N2) - (x2 * ones(1,N1))';
    x_log = - x_d .^2 / t2;
    GK = t1 * ((exp(1) * ones(N1,N2)) .^ x_log);
end

function LLH = LogLikelihood(x,y,tt1,tt2,tt3)
    Nx = length(x);
    K_theta = GaussianKernelMat(x, x, tt1, tt2) + tt3 * eye(Nx);
%    LLH = - log(det(K_theta)) - detrend(y') * inv(K_theta) * detrend(y);
    LLH = - log(det(K_theta)) - y' * inv(K_theta) * y;
end

5. 参考文献

『ガウス過程と機械学習』持橋大地・大羽成征 MLP
- 本ページトップにリンクあり
Qiita記事『pythonでガウス過程を実装する。パラメータの調整はマルコフ連鎖モンテカルロ法（MCMC）を使う』
- 今回はこのphyblasさんという方とほぼおんなじことをやっています。この方はpythonで書いてらっしゃいます。 qiita.com
『日本統計学会創立 75 周年記念出版 21 世紀の統計科学 < Vol. III > 数理・計算の統計科学』国友直人・山本拓監修東京大学出版会
- ネットで無料ダウンロードできます。素晴らしいpdfです。11章のMCMCの説明を参照しました。

6. さいごに・感想

Gaussian processを実装してすこしでも気持ちをわかろうと努力してみました。回帰では、結局kernel函数の選択がモデル選択に相当すると思われるので、結構重要ですね。

今回はは全部 $\textbf{y}$ の平均を $\textbf{0}$ でやってましたが、これは要は回帰の精度に自信のないデータから遠く離れた点はなんでもかんでも平均ゼロを予測するよーということで、モデルによってはあまり正しくないといえますので、改善の余地があります。

次回はではハイパーパラメータの決定として勾配法でやってみようかと思います。またデータが増えるごとに精度が良くなる様子をみるため、アニメーション化も行います(予定)。

今回はJulia素人が頑張って作った感満載で、また数式やコードじたいも本格的にブログに書いたのでえっらい時間がかかりました...... 研究します......

Gotterdammerung

JuliaによるGauss過程回帰 (MCMCでハイパーパラメータ決定)