@i_kiwamu Bayes - -

Size: px

Start display at page:

Download "@i_kiwamu Bayes - -"

まれあねぎたや
7 years ago
Views:

1 Bayes RStan

2 @i_kiwamu Bayes - -

3 Stan / RStan Bayes Stan Development Team - Andrew Gelman, Bob Carpenter, Matt Hoffman, Ben Goodrich, Michael Malecki, Daniel Lee and Chad Scherrer Open source Stand alone Stan R RStan

4 Stan by Andrew Gelman 1. graphical model compiler 2. (Adaptive) Hamiltonian Monte Carlo sampling 3. Gibbs sampling 4. etc......? ( д )

5 Stan BUGS JAGS Hamiltonian Monte Carlo No-U-Turn sampling...

6 1. RStan 2. Gibbs sampler

7 R Stan JAGS r_res_meeting_2012/ Mercurial GitHub

8 RStan Stan R 1. Running Stan from R 2. RStan Getting Started Guide

9 R Windows R install.packages(c( inline, Rcpp, RcppEigen ))

10 C++ compiler - Mac: Xcode - Windows: GCC in Rtools - Linux: GCC or Clang C++ compiler - Makeconf - Before: CXXFLAGS = -g -O2 $(LTO) - After: CXXFLAGS = -g -O3 $(LTO)

11 rstangettingstarted/8schools.stan data { int<lower=0> J; // number of schools real y[j]; // estimated treatment effects real<lower=0> sigma[j]; // s.e. of effect estimates parameters { real mu; real<lower=0> tau; real eta[j]; transformed parameters { real theta[j]; for (j in 1:J) theta[j] <- mu + tau * eta[j]; model { eta ~ normal(0, 1); y ~ normal(theta, sigma); BUGS

12 rstangettingstarted/8schools.stan data { int<lower=0> J; // number of schools real y[j]; // estimated treatment effects real<lower=0> sigma[j]; // s.e. of effect estimates parameters { real mu; real<lower=0> tau; real eta[j]; transformed parameters { real theta[j]; for (j in 1:J) theta[j] <- mu + tau * eta[j]; model { eta ~ normal(0, 1); y ~ normal(theta, sigma); Stan <lower=**,upper=**>

13 rstangettingstarted/8schools.stan data { int<lower=0> J; // number of schools real y[j]; // estimated treatment effects real<lower=0> sigma[j]; // s.e. of effect estimates parameters { real mu; real<lower=0> tau; real eta[j]; transformed parameters { real theta[j]; for (j in 1:J) theta[j] <- mu + tau * eta[j]; model { eta ~ normal(0, 1); y ~ normal(theta, sigma);

14 rstangettingstarted/8schools.stan data { int<lower=0> J; // number of schools real y[j]; // estimated treatment effects real<lower=0> sigma[j]; // s.e. of effect estimates parameters { real mu; mu tau real<lower=0> tau; real eta[j]; transformed parameters { uniform real theta[j]; for (j in 1:J) theta[j] <- mu + tau * eta[j]; msg/stan-users/4gv3fncqsnk/ model { VonPUkdcZuUJ eta ~ normal(0, 1); y ~ normal(theta, sigma); stan-reference pdf section C.3

15 rstangettingstarted/8schools.stan data { int<lower=0> J; // number of schools real y[j]; // estimated treatment effects real<lower=0> sigma[j]; // s.e. of effect estimates parameters { real mu; real<lower=0> tau; real eta[j]; transformed parameters { real theta[j]; for (j in 1:J) theta[j] <- mu + tau * eta[j]; model { eta ~ normal(0, 1); y ~ normal(theta, sigma); for 1 { 2

16 rstangettingstarted/8schools.stan data { int<lower=0> J; // number of schools real y[j]; // estimated treatment effects real<lower=0> sigma[j]; // s.e. of effect estimates parameters { real mu; real<lower=0> tau; real eta[j]; transformed parameters { real theta[j]; for (j in 1:J) theta[j] <- mu + tau * eta[j]; model { eta ~ normal(0, 1); y ~ normal(theta, sigma); dnorm normal tau sigma for

17 R > library(rstan) > schools.dat <- list(j = 8, y = c(28, 8, -3, 7, -1, 1, 18, 12), sigma = c(15, 10, 16, 11, 9, 11, 10, 18)) > fit1 <- stan(file='8schools.stan', data=schools.dat, iter=10000, chains=4, verbose=false) > print(fit1) > plot(fit1) > # extract data > fit1.la <- extract(fit1, permuted=true) > fit1.coef$mu <- median(fit1.la$mu)

18 R rstan > library(rstan) > schools.dat <- list(j = 8, y = c(28, 8, -3, 7, -1, 1, 18, 12), sigma = c(15, 10, 16, 11, 9, 11, 10, 18)) > fit1 <- stan(file='8schools.stan', data=schools.dat, iter=10000, chains=4, verbose=false) > print(fit1) > plot(fit1) > # extract data > fit1.la <- extract(fit1, permuted=true) > fit1.coef$mu <- median(fit1.la$mu)

19 R RStan > library(rstan) > schools.dat <- list(j = 8, y = c(28, 8, -3, 7, -1, 1, 18, 12), sigma = c(15, 10, 16, 11, 9, 11, 10, 18)) > fit1 <- stan(file='8schools.stan', data=schools.dat, iter=10000, chains=4, verbose=false) > print(fit1) > plot(fit1) > # extract data > fit1.la <- extract(fit1, permuted=true) > fit1.coef$mu <- median(fit1.la$mu)

20 R stan stan > library(rstan) > schools.dat <- list(j = 8, y = c(28, 8, -3, 7, -1, 1, 18, 12), sigma = c(15, 10, 16, 11, 9, 11, 10, 18)) > fit1 <- stan(file='8schools.stan', data=schools.dat, iter=10000, chains=4, verbose=false) > print(fit1) > plot(fit1) > # extract data > fit1.la <- extract(fit1, permuted=true) > fit1.coef$mu <- median(fit1.la$mu)

21 R print > library(rstan) > schools.dat <- plot list(j = 8, y = c(28, 8, -3, 7, -1, 1, 18, 12), sigma = c(15, 10, 16, 11, 9, 11, 10, 18)) > fit1 <- stan(file='8schools.stan', data=schools.dat, iter=10000, chains=4, verbose=false) > print(fit1) > plot(fit1) > # extract data > fit1.la <- extract(fit1, permuted=true) > fit1.coef$mu <- median(fit1.la$mu)

22 R extract > library(rstan) > schools.dat <- list(j = 8, y = c(28, 8, -3, 7, -1, 1, 18, 12), sigma = c(15, 10, 16, 11, 9, 11, 10, 18)) > fit1 <- stan(file='8schools.stan', data=schools.dat, iter=10000, chains=4, verbose=false) > print(fit1) > plot(fit1) > # extract data > fit1.la <- extract(fit1, permuted=true) > fit1.coef$mu <- median(fit1.la$mu)

23 plot(fit1) Rhat log probability

24 rstan::traceplot(fit1) burn-in

25 Gibbs sampler Nile - dlm package (1.1-0) - JAGS (3.3.0) - RStan (1.0.3) MCMC - iteration: burn-in: thinning rate: 1 - chain: 4 - Mac OS X R

26 Yi θi F i θi θi-1 G Y i = F i + i i N(0, y) i = G i 1 + i i N(0, )

27 dlm package > library(dlm) > nile.dlm <- dlmgibbsdig(nile, mod=dlmmodpoly(1), shape.y=0.1, rate.y=0.1, shape.theta=0.1, rate.theta=0.1, n.sample=8000) > burn <- 1:4000 > sigma.dlm <- matrix(unlist(nile.dlm[c("dv", "dw")]), ncol=2)[-burn,] > sigma.dlm <- sqrt(sigma.dlm) > colnames(sigma.dlm) <- c("sigma y", "sigma theta") > sigma.mcmc.dlm <- mcmc.list(mcmc(sigma.dlm)) > plot(sigma.mcmc.dlm)

28 dlm package σy σθ dlm chain 1

29 JAGS model{ for(i in 2:N){ Nile[i] ~ dnorm(f * theta[i], tau.y) theta[i] ~ dnorm(g * theta[i-1], tau.theta) Nile[1] ~ dnorm(f * theta[1], tau.y) theta[1] ~ dnorm(nile.1, tau.theta) F ~ dnorm(0,.0001) G ~ dnorm(0,.0001) tau.y <- pow(sigma.y, -2) tau.theta <- pow(sigma.theta, -2) sigma.y ~ dunif(0, 1000) sigma.theta ~ dunif(0, 1000)

30 JAGS > library(r2jags) > nile.1 <- Nile[1] > nile.data.jags <- list("n", "Nile", "nile.1") > nile.inits.jags <- function() { list(theta=rnorm(n), sigma.theta=runif(1), sigma.y=runif(1), F=rnorm(1), G=rnorm(1)) > nile.params.jags <- c("theta", "sigma.y", "sigma.theta", "F", "G") > nile.jags <- jags(nile.data.jags, nile.inits.jags, nile.params.jags, model.file="nile.jags", n.iter=8000, n.burnin=4000, n.chains=4) > sigma.mcmc.jags <- lapply(1:4, function(i) as.mcmc(nile.jags)[[i]][,c(5,4)]) > plot(sigma.mcmc.jags)

31 JAGS

32 RStan data { int<lower=0> N; real<lower=0> Nile[N]; parameters { real theta[n]; real F; real G; real<lower=0> sigma_y; real<lower=0> sigma_theta; model { theta[1] ~ normal(nile[1], sigma_theta); Nile[1] ~ normal(f * theta[1], sigma_y); for(i in 2:N){ theta[i] ~ normal(g * theta[i-1], sigma_theta); Nile[i] ~ normal(f * theta[i], sigma_y);

33 RStan > nile.data.stan <- list(n=n, Nile=Nile) > nile.stan <- stan(file="nile.stan", data=nile.data.stan, iter=8000) > nile.stan.extract <- extract(nile.stan, inc_warmup=false) > sigma.mcmc.stan <- mcmc.list(mcmc(nile.stan.extract[,1,103:104]), mcmc(nile.stan.extract[,2,103:104]), mcmc(nile.stan.extract[,3,103:104]), mcmc(nile.stan.extract[,4,103:104])) > plot(sigma.mcmc.stan)

34 RStan

35 RStan observed estimation 95% HPD

36 RStan σy 95% σθ 95% Iteration 1 50 chain 1: chain 2: chain 3: chain 4:

37 RStan σy 95% Iteration chain 1 3 chain 4 σθ 95% chain 1: chain 2: chain 3: chain 4:

38 RStanの収束の様子 σの 95%範囲収束完了 200 Iteration chain 1: 青 σθの 95%範囲 0 sigma theta y sigma y chain 2: 緑 chain 3: 黄 chain 4: 赤

39 RStan JAGS θ3の初期値 Burn-in期間中に 500 θ3 JAGS RStanで再計算 1500 同じ初期値を設定し 1000 幅広いサンプリング推定結果が初期値に左右されにくい 0 幅広くサンプリング θ2の初期値 θ

40 150 / dlm package JAGS RStan RStan JAGS 3

41 data { int<lower=0> N; real y[n]; transformed data { real sum_log_y; sum_log_y <- 0.0; for(i in 1:N) sum_log_y <- sum_log_y + log(y[i]); parameters { real<lower=0> mu; real<lower=0> lambda; f (x µ, ) = 2 x 3 model { lp <- lp * N * log(lambda / (2.0 * pi())) * sum_log_y; for(i in 1:N) lp <- lp - lambda * pow((y[i] - mu) / mu, 2) / (2.0 * y[i]); 1/2 exp (x µ) 2 2µ 2 x lp

42 > library(suppdists) > y <- rinvgauss(100, nu=5, lambda=50) > inv.gauss.data <- list(n = 100, y = y) > inv.gauss.stan <- stan(file = inv_gauss.stan, data = inv.gauss.data) > print(inv.gauss.stan) Inference for Stan model: inv_gauss. 4 chains: each with iter=2000; warmup=1000; thin=1; 2000 iterations saved. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat mu lambda lp Samples were drawn using NUTS2 at Sun Nov 25 17:04: For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).

43 Stan/RStan Gibbs sampler - - -

44 Stan reference stan-reference pdf ML or - ML: hl=ja&fromgroups#!forum/stan-users - FAQ: stan-users/4gv3fncqsnk/x95dtubcgcwj

45 Thank you for your attentions!

2014ESJ.key

2014ESJ.key http://www001.upp.so-net.ne.jp/ito-hi/stat/2014esj/ Statistical Software for State Space Models Commandeur et al. (2011) Journal of Statistical Software 41(1) State Space Models in R Petris & Petrone (2011)