2014-03-28

rstan

BUGS ベイズ MCMC R stan rstan Gibbs sampling

Stanについてはこちらに紹介記事があり、環境整備に関しては、こちらを参考にして、以下のようにやればよいようです(Windowsの場合)
- 基本的には
  - Rがあって
  - Rtoolsを入れて、Rcppパッケージが動くようにして、inlineでのRcppによるRのC++化を可能にし
  - rstanを特別なところからダウンロードして、という流れです。
- 以下、その方法ソース

# Rtoolsをインストール＠https://github.com/stan-dev/rstan/wiki/Install-Rtools-for-Windows
Sys.getenv('PATH')

# Rcppが動くかどうかの確かめ
install.packages('inline')
install.packages('Rcpp')
# using library inline to compile a C++ code on the fly
library(inline) 
library(Rcpp)
src <- ' 
  std::vector<std::string> s; 
  s.push_back("hello");
  s.push_back("world");
  return Rcpp::wrap(s);
'
hellofun <- cxxfunction(body = src, includes = '', plugin = 'Rcpp', verbose = FALSE)
cat(hellofun(), '\n') 

# rstanをインストール。ただしいわゆるCRANからのインストールではないので、その指定をしつつ…
Sys.setenv(R_MAKEVARS_USER = "foobar")
options(repos = c(getOption("repos"), rstan = "http://wiki.rstan-repo.googlecode.com/git/"))
install.packages('rstan', type = 'source')
Sys.unsetenv("R_MAKEVARS_USER")

# rstanを使うための２つの条件設定
library(rstan)
set_cppo("fast")  # for best running speed

library(rstan)
set_cppo('debug') # make debug easier

その上でこちらをなぞってみましたが、コンパイルエラー…。このあたりに弱い私は、ここで面倒臭くなりました
- 将来、戻ってくるかもしれないので、そのソースとファイルを以下に張り付けておきます

schools_code <- '
  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);
  }
'

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))

fit <- stan(model_code = schools_code, data = schools_dat, 
            iter = 1000, chains = 4)
            
fit1 <- stan(file = '8schools.stan', data = schools_dat, 
             iter = 1000, chains = 4)
fit2 <- stan(fit = fit1, data = schools_dat, iter = 10000, chains = 4)

print(fit2)
plot(fit2)

la <- extract(fit2, permuted = TRUE) # return a list of arrays 
mu <- la$mu 

### return an array of three dimensions: iterations, chains, parameters 
a <- extract(fit2, permuted = FALSE) 

### use S3 functions as.array (or as.matrix) on stanfit objects
a2 <- as.array(fit2)
m <- as.matrix(fit2)

- 呼出しstanファイル "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);
}
||


*1395781231*[BUGS][ベイズ][MCMC][openBUGS][Gibbs sampling]openBUGS
-Bayesian inference using Gibbs sampling:BUGS は、ベイズ推定を計算機統計学的にやってしまいましょうという手法
-計算機統計学的に…という部分は、「難しい関数記述(の省略できるところ)省略して、MCMC/Gibbs samplingで代用する」という風に読み替えてもよい
-Rと連動する使い方として、DynacomさんのサイトにBUGSを使う場合([http://www.dynacom.co.jp/tech/tech_column/tech009.html:title=こちら])とStanを使う場合([http://www.dynacom.co.jp/tech/tech_column/tech013.html:title=こちら])との紹介があります
*1395781232*[BUGS][ベイズ][MCMC][R][stan][rstan][Gibbs sampling]rstan
-Stanについては[http://tjo.hatenablog.com/entry/2014/03/25/194605:title=こちら]に紹介記事があり、環境整備に関しては、こちらを参考にして、以下のようにやればよいようです(Windowsの場合)
--基本的には
---Rがあって
---Rtoolsを入れて、Rcppパッケージが動くようにして、inlineでのRcppによるRのC++化を可能にし
---rstanを特別なところからダウンロードして、という流れです。
--以下、その方法ソース
>|r|
# Rtoolsをインストール＠https://github.com/stan-dev/rstan/wiki/Install-Rtools-for-Windows
Sys.getenv('PATH')

# Rcppが動くかどうかの確かめ
install.packages('inline')
install.packages('Rcpp')
# using library inline to compile a C++ code on the fly
library(inline) 
library(Rcpp)
src <- ' 
  std::vector<std::string> s; 
  s.push_back("hello");
  s.push_back("world");
  return Rcpp::wrap(s);
'
hellofun <- cxxfunction(body = src, includes = '', plugin = 'Rcpp', verbose = FALSE)
cat(hellofun(), '\n') 

# rstanをインストール。ただしいわゆるCRANからのインストールではないので、その指定をしつつ…
Sys.setenv(R_MAKEVARS_USER = "foobar")
options(repos = c(getOption("repos"), rstan = "http://wiki.rstan-repo.googlecode.com/git/"))
install.packages('rstan', type = 'source')
Sys.unsetenv("R_MAKEVARS_USER")

# rstanを使うための２つの条件設定
library(rstan)
set_cppo("fast")  # for best running speed

library(rstan)
set_cppo('debug') # make debug easier

その上でこちらをなぞってみましたが、コンパイルエラー…。このあたりに弱い私は、ここで面倒臭くなりました
- 将来、戻ってくるかもしれないので、そのソースとファイルを以下に張り付けておきます

schools_code <- '
  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);
  }
'

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))

fit <- stan(model_code = schools_code, data = schools_dat, 
            iter = 1000, chains = 4)
            
fit1 <- stan(file = '8schools.stan', data = schools_dat, 
             iter = 1000, chains = 4)
fit2 <- stan(fit = fit1, data = schools_dat, iter = 10000, chains = 4)

print(fit2)
plot(fit2)

la <- extract(fit2, permuted = TRUE) # return a list of arrays 
mu <- la$mu 

### return an array of three dimensions: iterations, chains, parameters 
a <- extract(fit2, permuted = FALSE) 

### use S3 functions as.array (or as.matrix) on stanfit objects
a2 <- as.array(fit2)
m <- as.matrix(fit2)

- 呼出しstanファイル "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);
}