2020-03-07

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

初っ端ですが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素人が頑張って作った感満載で、また数式やコードじたいも本格的にブログに書いたのでえっらい時間がかかりました...... 研究します......

2020-02-20

Blogを再開しようと思う

Blogを開設すると言って、2回更新しただけで、早くも3年(笑)が経ってしまいました。

私はいまどこで何をしているかというと、沖縄県の病院から東京の大学に戻り、修士1年の院生をやり直しています。
東京大学の総合文化研究科というところで、医学系の研究科ではありません。
研究内容は理論神経科学、ラボのメインテーマは「意識」です。

なぜ臨床医をやめて修士に入学し直したのか？研究といっても医学系でないのはなぜか？については別postで述べようと思います。
意識については中学生くらいから興味を持っていました。
医師をやりながら、いつかは意識の研究に取り組みたいとは前々から思っていたのですが、それが意外と早くこの時期だったと言うことです。

Blogを再開するのは次の理由からです;
① Outputの機会を作る
② 自分に課す
③ 備忘録

① Outputの機会を作る

学会発表やインタビューや研究期の応募文など、アウトプットの機会はこれから無数にあるわけですが、思えばきちんとアウトプットの練習をする機会ってないなぁと思うわけです。
思えば自分もまったくアウトプットが得意ではなく、高校時代の国語の成績などは壊滅的でした。
仕事ができればなんでもいいだろうと思った時期もありましたが、自分の考えをまとめて発表するということは練習を積まないと出来ないし、
優れた研究者や友人でも尊敬できる方は積極的にアウトプットをしています。
ということで、自分もやってみようと思います。

② 自分に課す

自分のようなぐうたら人間は目標や締め切りが決まっていないと何事も腰が重い。
途中まで読んでほっぽりだした教科書など数知れず……
現状としても、大学院に入って半年、情報系や応用数学の知識や技術は以前よりは身についたとは思いますが、自分の中でのブレイクスルーを見いだしていません。
そこで、学習記録として自分にこのブログを課すことで、自分に葉っぱをかけようじゃあないかと言うわけです。
内容は現在は統計・機械学習の勉強・実装記録をメインにしようと思っていますが、教科書や論文のまとめ、読書感想文、考えたことなどもポストする予定です。
とりあえず、Netflix「テラスハウス」と同じ頻度、(最低)月3回を目指します。

③ 備忘録

あまりこれまで記録というものをつけてこなかった人生で、単純に勉強したことや考えたことの記録を作りたいというのもあります。

よろしくお願いします。

2017-04-25

学生時代の勉強について各論

アメリカを目指す全ての医学生が立ち向かわなければいけない壁、USMLE.

先日3/27にStep 1を受験し、261点で合格しました。

誰かの役に立てば良いと思い、勉強に関して考えることをいまのうちに書いておこうと思います。

質問ある方は大歓迎です。

勉強総論と勉強各論を考えたのですが、総論を書く前に各論がかけてしまったのであげておきます笑

学生時代の医学の勉強に関して、経験から次のことが言えると思います。

☆ 勉強の3要素

① 新知識のインプット

② 既習知識のアウトプット

③ 忘却を最小限に抑えること

そしてこの3つそれぞれに対応する教材が存在します。僕はこの3つを中核として勉強を進めました。

① => The First Aid for USMLE Step 1 (2015, 2016, 2017)

② => Online qbanks (USMLE world, Kaplan qbank)

③ => Anki

以下詳述します。

① The First Aid for USMLE Step 1

USMLE Step 1を受験するならばなしには語れない本。かつ、筆者の現在のあらゆる医学的知識の礎となっている本。これほど簡潔かつ網羅的に疾患の病態生理、薬の作用機序を説明した本を私は知りません。まさにbeauty.

ただ、最も大切な病態のメカニズムが本当にさらっと、ひとことで片付けられていたりするので、初学者は読むのに苦労するかもしれません。Step 1を受験しない人も、前半の基礎医学について読むと内科学の理解の度合いがぐんと上がると思います。個人的にmicrobiologyの章は白眉。

Webでよく見る勉強の王道は、このFAの余白に知識を書き込むこと、そしてそれを何度も読み、内容を全て覚えてしまうと言うことです。僕も結局通しで3-4周読んだと思いますが、後述するAnkiを使ったこともあり、この本に書かれている内容は結構覚えていたと思います。

しかし、僕個人的には書き込み→読み込むという勉強法はあまり合っていなかった。これは人によると思うんですが、僕は新しい知識を勉強を本に書き込むと、なんだかそれだけで満足してしまいあまり読み返さない人でした。あと、自分の手書きの文字が内容としてパッと入って来ず、活字より読みにくいのであまり意味がなかったのです。

加えて、この本に頼りすぎない、この本に振り回されないのも重要です。5年次の受験予定日直前にはこの本をはしからはしまで暗記しようと思って読んでいましたが、模試の点数はあがりませんでした。重箱の隅をつつくようなfactもたまに書いてありますが、はっきり言って殆ど出ません。

それよりも基本的な事項を理由付けして、病態生理をきちんとおって、それごと覚えてしまう、という方が良いでしょう。闇雲に暗記するのではなく、なぜか？を常々考えましょう。

eg. Alzheimer disease では τ proteinの沈着が見られる。

これを記憶して見ましょう。

ダメな例)

Alzheimer ではτ protienか。そーゆーものなんだ。暗記暗記。

良い例)

Τ proteinの役割は細胞内微小管の安定化。Alzheimerでは神経細胞の中のτ proteinのhyperphosphorylationが問題となる。過リン酸化されるとτ protein は paired helical filament (PHF) となり、これがAlzheimerの神経細胞に見られるneurofibrially tangles の原因となる。PHFとなったτ protein は微小管をうまく安定化できず、結果神経細胞の細胞内輸送が障害される。だから神経伝達物質の輸送がうまくいかず認知症になる！！

この例をみてもわかるように、きっちり理由まで含めて覚えるのははっきり言って面倒くさいです。でも少しここで遠回りをしておけば、必要事項を覚えやすく、かつ思い出しやすくなりますし、似たような事項とも関連づけるのが容易です(この例だとたとえば、frontotemporal dementiaは同じような機序でおこるので、τ proteinが見える、とか)。

② Online q-banks

USMLEのような試験では、①のinputだけ増やしてもどうしようもないところがあり、実際に問われるのは知識の多さではなく、知識の運用法である。知識のアウトプットの練習として、問題をたくさん解くことは非常に(ある意味最も)重要。

Q-bankとしては、USMLE World, Kaplan qbanks, USMLERxの3つが有名。難易度的にはUW ≧ Kaplan >> USMLERx といったところか。あと、USMLEの作問機関が作っているNBME, UWの作っているUWSAという模試もやると力試しになります。

どれに取り組むにせよ、解説まできちんと読み、重要事項や気づかなかった点をAnkiカードにしてしまうという勉強法が良い気がする(Ankiカードの作り方についてはAnkiの項で詳述)。

5年次ではUSMLERx 1周、USMLE World 3周、Kaplan qbanks 0.8周に取り組みました。今思うとUWを3周したのは、問題ごと覚えていた問題も多くあまり意味がなかったとは思います。2周くらいでいいんじゃないでしょうか。NBMEの模試は4つくらい受け、得点は

NBME11: 228

NBME15: 243

NBME13: 237

NBME17: 237

でした。

6年次にはUSMLE Worldをもう1周、Kaplanを1周しました。模試は3つ受け、

UWSA1: 264

NBME18: 246

UWSA2: 251

でした。

③ Anki

暗記補助アプリ。簡単にいうとiphone版フラッシュカードで、フラッシュカードの問題に対して、覚えていないカード、答えられるが身についていないカード、普通に覚えたカード、楽勝なカード、から選択することで、カードの学習間隔が変わっていきます。この方法であらゆる知識を長期記憶として脳に保存し、「忘れることを、忘れる」ことができるわけです。

具体的には、FAや自分で勉強したfactsをカードとして加えていけばよし。コツとしては、①でも述べましたが、なるべく1問1答のようなカードを作らないこと。病態生理まで暗記できるようなカードを作ること。例えばさっきの①の例だと、

ダメなAnkiカードの例)

———————————————————————————————————

(表)

What can be seen in Alzheimer disease microscopically?

———————————————————————————————————

(裏)

Tau protein

———————————————————————————————————

良いAnkiカードの例)

———————————————————————————————————

(表)

What is the function of τ protein? Associated diseases (2)

Provided this function, what can be explained?

———————————————————————————————————

(裏)

Tau protein: binds to and stabilizes microtubles

Alzheimer dementia (AD) & frontotemporal dementia

If tau proteins are hyperphosphorylated -> paired helical filament (PHF) tau = major component of neurofibrially tangles -> CANNOT bind to microtubules so impairment in neuronal transport systems -> # of tangles associated /c degree of disease

UpToDate:

Tau is a microtubule-associated protein that aids in microtubule assembly and stabilization.

In AD, tau becomes hyperphosphorylated and aggregates to form paired helical filament (PHF) tau, a major component of neurofibrillary tangles within the neuronal cytoplasm. The accumulation of this altered protein is toxic to neurons in experimental models. In addition, transmission of pathologic forms of tau between neurons has been proposed to account for the spread of AD in brain, which follows a distinct progression across brain regions as AD advances

———————————————————————————————————

良い例では、

(1) 病態生理に重きを置き、理由付けをしている

(2) 改行、太字などを用いて覚えるべきところをわかりやすくしている

ことがわかりますね。この方が断然覚えやすいです。

こんなふうにしてqbank, FAの内容などをカードにして覚えていきましょう。

また、Ankiに関してはYousmle.comのAlecの記事が非常に参考になります。

その他

そのほかに手をつけた本も上げて置きます。ほとんどは通読していない、つまみ食いだけです。

● Rapid Review of Pathology (Goljan)

American 三苫, American KSR, 呼び名はいろいろありますが、米国のUSMLE のカリスマ講師、ゴルジャン先生直筆の教科書。FAの副読本的立ち位置と勝手に思っている。結構細かく病態生理まで書いてあって参考になる。箇条書きなので知識が繋がりにくいのが個人的には難。

● Biochemistry (Lippincotr’s Illustrated Reviews)

生化学の読みやすい教科書。生化学ははじめとっつきにくいですが、訳がわかってしまうと楽しくなります。僕はFAの生化学の項と、必要に準じてこの本を参照していました。図がわかりやすいので図ごとAnkiにつっこんでしまってもいいかもしれない。

● Principles of Pharmacology (Golan)

ハーバードで薬理の教科書として使われているという例の本。たしかに学生で必要な薬理の知識が詳しく載っていてためになる。図もわかりやすい。

● Physiology (Linda Costanzo)

生理学のイチオシ。医学部3年のときにこの日本語版を読んでいましたが、説明がクリアカットでとてもわかりやすい。初学者には最適。ただクリアカットすぎて、たまに重要事項が抜けていたりする。けれどはじめは問題にならないので良いと思います。

Eg. VSDではcyanosisにならない。=> 実際はなる。欠損孔が大きいとLVから肺に駆出される血液が多くなりV/Q mismatchを引き起こすため。

● USMLE Step 2 CK Secrets

国試・マッチング対策を日本語でするのがいやで読んでいた本。Step 1に関しても、最近臨床的な問題が増えているから、Step 2 CKの知識が役に立つ部分は少なくないと思います。

● 研修医当直御法度第6版

通称赤本。研修医が救急科当直で出会いそうなsituationの対応と陥りがちなpitfallについて記してある。とてもわかりやすく書かれており、1ページあたりの内容も多くはないので、一冊丸ごと覚えるような勢いで読んでもいいかも知れない。ここで読んだ知識は国家試験の外科・救急領域でもかなり役に立ちました。

● マッチングの過去問

僕の志望した病院は、全問英語でしかも試験の内容も難しかったので、過去問を勉強するだけでとても勉強になりました。僕の頃は筆記+面接だったのですが今年から面接だけになってしまうようで少し心配です。

以下、読んでないけど参考になるかも知れない本

● 内科レジデントの鉄則……聖路加の教育法

● OCH救急マニュアル……沖縄中部の教育法

● Step beyond Resident……読み物で読みやすい

● キクタン……日本語で医学的知識がある程度ある人におすすめ

ご意見などありましたらお聞かせください。

2017-04-19

はじめまして

みなさまはじめまして。

東京の医学部卒業後、2017年より沖縄県で初期研修医として働いているgvskです。

将来の進路として、米国での外科研修を経て、向こうで心臓血管外科医になりたいと考えています。

アメリカでの臨床研修に興味のある後輩の方々に少しでも参考になるかなと思ってブログを始めました。

医師になりたての僕ですが、これから自分がどうなるのか自分でもわからないので、その備忘録も兼ねています。

経歴:

神奈川県S高校卒業

↓

東京の大学の医学部卒業

↓

沖縄の病院で初期研修

趣味はピアノ、テニス、オペラ鑑賞(特にWagner)、料理、読書、車の運転、筋トレ

ブログには趣味の話も載せられればと思います。

需要があるかわかりませんが、なにとぞ宜しくお願いします。