ベイズ塾の春の合宿で発表した資料です。
スライドシェアにアップしました。
ガウス過程についての解説と、Stanでの実行方法などについて書きました。
追記：ガウス過程のコードは松浦さんのサイトを参考にしています。ありがとうございます。
※一部数式が間違えてたので修正しました。
スライドで使っているデータとコードを以下にアップします。
データ
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
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コード
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 |
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 |