Rでガウス過程回帰の基礎

# <font color="DarkSlateBlue">記事の目的</font>

この記事では、ガウス過程でいろいろなカーネルを使ってサンプリングします。
また、予測分布も可視化してみます。
参考: [ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)](https://www.amazon.co.jp/%E3%82%AC%E3%82%A6%E3%82%B9%E9%81%8E%E7%A8%8B%E3%81%A8%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92-%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%83%97%E3%83%AD%E3%83%95%E3%82%A7%E3%83%83%E3%82%B7%E3%83%A7%E3%83%8A%E3%83%AB%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA-%E6%8C%81%E6%A9%8B-%E5%A4%A7%E5%9C%B0/dp/4061529269)


# <font color="DarkSlateBlue">目次</font>

<table>
<tr><td><font color="DarkslateBlue">No.<td><font color="DarkslateBlue">目次
<tr><td><font color="DarkslateBlue">1</font><td><a href="#1-いろいろなカーネルの計算">いろいろなカーネルの計算</a>
<tr><td><font color="DarkslateBlue">2<td><a href="#2-ガウス過程の可視化">ガウス過程の可視化</a>
<tr><td><font color="DarkslateBlue">3<td><a href="#3-ガウス過程回帰の予測分布">ガウス過程回帰の予測分布</a>
<tr><td><font color="DarkslateBlue">4<td><a href="#4-予測分布の可視化">予測分布の可視化</a>
</table>

# <font color="LimeGreen">1.</font> <font color="DarkSlateBlue">いろいろなカーネルの計算</font>
![IMG_0202.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/775500/88e00ae9-688a-0668-11c0-b48be4cfa711.jpeg)



```R:
#データ
X <- seq(-4,4, 0.1)
N <- length(X)
#ガウスカーネル計算
K1 <- matrix(NA, nrow=N, ncol=N)
K2 <- matrix(NA, nrow=N, ncol=N)
K3 <- matrix(NA, nrow=N, ncol=N)
K4 <- matrix(NA, nrow=N, ncol=N)
theta1.1 <- 1
theta1.2 <- 1
theta3 <- 1
theta4.1 <- 1
theta4.2 <- 0.5
for(i in 1:N){
  for(j in 1:N){
    #ガウスカーネル計算
    K1[i,j] <- exp(-(X[i]-X[j])^2/theta1.2)/theta1.1
    #線形カーネル計算 (1次元ならfor文必要なし)
    K2 <- X%*%t(X)
    #指数カーネル計算
    K3[i,j] <- exp(-abs(X[i]-X[j])/theta3)
    #周期カーネル計算
    K4[i,j] <- exp(theta4.1*cos(abs(X[i]-X[j])/theta4.2))
  }
}
```

# <font color="LimeGreen">2.</font> <font color="DarkSlateBlue">ガウス過程の可視化</font>
```R:
library(mvtnorm)
library(ggplot2)
library(gridExtra)
set.seed(10)
kernel.plot <- function(K, name){
  y1 <- rmvnorm(1, rep(0,N), K)
  y2 <- rmvnorm(1, rep(0,N), K)
  y3 <- rmvnorm(1, rep(0,N), K)
  y4 <- rmvnorm(1, rep(0,N), K)
  ggplot(NULL,aes(x=X)) +
    geom_line(aes(y=y1)) + geom_line(aes(y=y2)) + 
    geom_line(aes(y=y3)) + geom_line(aes(y=y4)) + 
    ylim(-6, 6) + labs(x="x",y="y", title=name)
}
p1 <- kernel.plot(K1, "ガウスカーネル")
p2 <- kernel.plot(K2, "線形カーネル")
p3 <- kernel.plot(K3, "指数カーネル")
p4 <- kernel.plot(K4, "周期関数")
grid.arrange(p1, p2, p3, p4)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/775500/d1f0f84f-197f-d810-54b9-aa1a18ec534b.png)

# <font color="LimeGreen">3.</font> <font color="DarkSlateBlue">ガウス過程回帰の予測分布</font>
![IMG_0201.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/775500/f64e54e2-7175-9ed6-f184-781186e61ac0.jpeg)


```R:
#データ
set.seed(100)
N <- 10
M <- 101
X <- append(rnorm(N, 0, 5), seq(-5,5,0.1))
#カーネルの計算
K <- matrix(NA, nrow=N+M, ncol=N+M)
theta1 <- 1
theta2 <- 1
for(i in 1:(N+M)){
  for(j in 1:(N+M)){
    #ガウスカーネル
    K[i,j] <- exp(-(X[i]-X[j])^2/theta2)/theta1
  }
}
K <- K + 0.1*diag(nrow(K))
y <- rmvnorm(1, rep(0,N+M), K)
#予測分布の計算
mu <- K[(N+1):(N+M),1:(N)] %*% solve(K[1:N,1:N]) %*% y[1:N]
var <- K[(N+1):(N+M), (N+1):(N+M)] - 
  K[(N+1):(N+M),1:(N)] %*% solve(K[1:N,1:N]) %*% K[1:N, (N+1):(N+M)]
```

# <font color="LimeGreen">4.</font> <font color="DarkSlateBlue">予測分布の可視化</font>
```R:
sigma.2 <- diag(var)
ggplot(NULL,aes(x=X[(1+N):(N+M)])) +
  geom_line(aes(y=mu), size=1.2, color="blue") + 
  geom_ribbon(aes(ymin=mu+2*sigma.2,ymax=mu-2*sigma.2), alpha=0.1, fill="blue") + 
  geom_point(aes(x=X[1:N],y=y[1:N]), color="blue", shape=1, size=5) +
  xlim(-5, 5)+ ylim(-2.5,2.5) +
  labs(x="x",y="y", title="ガウス過程回帰の予測分布")
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/775500/502ad603-7cd2-66bd-db0b-4fa6618732a1.png)