So-net無料ブログ作成
  • ブログをはじめる
  • ログイン

Stan: ガウス過程をためしてみたメモ [統計]

Stanでガウス過程をためしてみたメモです。

だいたいStan Modeling Language User’s Guide and Reference Manualの18章を参考にしました。また、 StatModeling Memorandumに、Stanによるガウス過程の記事があります。

Stanのコードは、Stan example-modelsgp-predict.stanを使用しました。

Rコードは以下のとおりです。

library(ggplot2)
library(dplyr)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

set.seed(108)
N <- 40
sd <- 0.2
x <- runif(N, 0, 2 * pi)
y <- sin(x) + rnorm(N, 0, sd)
df <- data.frame(x = x, y = y)
x_new <- seq(0.25 * pi, 1.75 * pi, length = 7)

model_file <- "https://raw.githubusercontent.com/stan-dev/example-models/master/misc/gaussian-process/gp-predict.stan"

fit <- stan(file = model_file,
            data = list(N1 = N,
                        x1 = x,
                        y1 = y,
                        N2 = length(x_new),
                        x2 = x_new))
print(fit, pars = c("rho", "alpha", "sigma", "lp__"))

y2_summary <- extract(fit, pars = "y2")[[1]] %>%
  apply(2, quantile, probs = c(0.025, 0.5, 0.975))
y2 <- data.frame(x = x_new,
                 low = y2_summary[1, ],
                 med = y2_summary[2, ],
                 high = y2_summary[3, ])

ggplot(y2) +
  geom_point(aes(x, med), colour = "red") +
  geom_linerange(aes(x = x, ymin = low, ymax = high),
                 colour = "red") +
  geom_point(data = df, aes(x, y)) +
  stat_function(fun = sin, linetype = 2)

結果です。

Rplot01.png


タグ:STAn
nice!(1)  コメント(0) 
共通テーマ:日記・雑感

nice! 1

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント