Stanでガウス過程

ベイズ塾の春の合宿で発表した資料です。

スライドシェアにアップしました。

ガウス過程についての解説と、Stanでの実行方法などについて書きました。

追記:ガウス過程のコードは松浦さんのサイトを参考にしています。ありがとうございます。

※一部数式が間違えてたので修正しました。

スライドで使っているデータとコードを以下にアップします。

データ

    x          y
1   1 0.07051282
2   2 0.12820513
3   3 0.13461539
4   4 0.16666667
5   5 0.21153846
6   6 0.23717949
7   7 0.24358974
8   8 0.26282051
9   9 0.26282051
10 10 0.28846154
11 11 0.30128205
12 12 0.31410256
13 13 0.33333333
14 14 0.33974359
15 15 0.35897436
16 16 0.39743590
17 17 0.40384615
18 18 0.42948718
19 19 0.44230769
20 20 0.48717949

 

Rコード

dat <- read.csv("gp_test.csv")

N1 <- nrow(dat)
x1 <- dat$x
y1 <- dat$y
N2 <- 59
x2 <- seq(1, 30, length.out = N2)

data.gp <- list(
  N1 = N1,
  N2 = N2,
  x1 = x1,
  x2 = x2,
  y1 = y1
)

model.gp <- stan_model("GP_reg.stan")

fit.gp <- sampling(
  model.gp,
  data = data.gp,
  iter = 5000,
  chains = 4,
  cores = 4
)

print(fit.gp, pars = "y2")
print(fit.gp, pars = "theta")

# plot
ext.fit <- rstan::extract(fit.gp)

y2.med <- apply(ext.fit$y2, 2, median)
y2.low <- apply(ext.fit$y2, 2, quantile, probs = 0.05)
y2.high <- apply(ext.fit$y2, 2, quantile, probs = 0.95)

d.est <-
  data.frame(x = x2,
             y = y2.med,
             ymax = y2.high,
             ymin = y2.low)
d.obs <-
  data.frame(
    x = x1,
    y = y1,
    ymax = rep(0, N1),
    ymin = rep(0, N1)
  )
p <- ggplot(d.est, aes(
  x = x,
  y = y,
  ymax = ymax,
  ymin = ymin
))
p <- p + xlab("x") + ylab("y")
p <- p + geom_ribbon(alpha = 1 / 3)
p <- p + geom_line(aes(y = y), color = "blue")
p <- p + geom_point(data = d.obs)
p

 

 

This entry was posted in 心理統計学. Bookmark the permalink.