2016-06-12

実行：ぱらぱらめくる『統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―』

実行は、srcAディレクトリのファイルをコピーペーストして実行するだけ
６章分

source('myfunc/myfunc.R')

冒頭でmyfuncディレクトリのmyfunc.Rソースを読み込みます
このmyfunc.Rは、myfuncディレクトリ内のそのほかのいろいろなxxx.Rソースを読み込みます
第１章用のファイル"part1.R"はstanを使うより前
第２章用から、stanを使います。それは第２章用ファイル"part2.R"の冒頭に

library(rstan)

としていることからわかります
どこでstanをしているかと言えば、part2.Rの中で使われている

#正規分布に関する推測
out <-G1mean(x,prior=T,mL=0, mH=1000, sL=0, sH=100)

にあるG1mean()という関数が見慣れません
これはmyfuncディレクトリ内にある関数で、以下のように中でrstanパッケージを呼び出したりしており、stanによる推定処理を決めています

##########################################################################
#1群の正規分布に関する推測
#■入力
#x:ベクトル形式のデータ 
#prior=F : 論理値。Tなら事前分布の範囲を指定。Fなら指定せず。
#mL, mH, sL, sH : 事前分布のパラメタ
#prob=c(0.025, 0.05, 0.5, 0.95, 0.975): 事後分布で報告する確率点
#see=1234,cha=5,war=1000,ite=21000 : MCMC関連のパラメタ
#■出力
#fit:stanの出力, par:母数リスト, prob:確率ベクトル, 
#mu:平均, sigma:標準偏差, xaste:予測分布, log_lik:対数尤度
##########################################################################
G1mean <- function(x,prior=F, mL=-1000, mH=1000, sL=0, sH=100,
prob=c(0.025, 0.05, 0.5, 0.95, 0.975),see=1234,cha=5,war=1000,ite=21000)
{
   library(rstan)
   if (prior) {
     dat <- list(n=length(x),x=x,mL=mL,mH=mH,sL=sL,sH=sH)
     scr <- "stan/G1meanPT.stan"                           # Stan ファイル名
   } else {
     dat <- list(n=length(x),x=x)
     scr <- "stan/G1meanPF.stan"                           # Stan ファイル名
   }
   par<-c("mu","sigma","xaste","log_lik")        # sampling変数
   fit<-stan(file=scr,data=dat,pars=par,seed=see,chains=cha,warmup=war,iter=ite)
   ext<-extract(fit, c("mu","sigma","xaste","log_lik"))
   mu<-ext$mu; sigma<-ext$sigma; xaste<-ext$xaste; log_lik<-ext$log_lik
   out<-list(fit=fit,par=par,prob=prob,mu=mu,sigma=sigma,xaste=xaste,log_lik=log_lik)
   class(out)<-'G1mean'
   return(invisible(out))
}

その内部で

scr <- "stan/G1meanPT.stan"                           # Stan ファイル名

という行がありますが、stanディレクトリ内の"G1meanPT.stan"というファイル(Stanのやり方を指定するファイル)を読み込んでいます
この.stanファイルは、以下のようになっていて、ある書式に従って書かれているらしきことが読み取れます。

data {
int<lower=0> n;                                  #データ数
real<lower=0> x[n];                              #データ
real mL; real mH; real sL; real sH;  #事前分布
}
parameters {
real<lower=mL,upper=mH> mu;             #平均(範囲指定)
real<lower=sL,upper=sH> sigma;          #標準偏差(範囲指定)
}
transformed parameters {
}
model {
x ~ normal(mu,sigma);                   #正規分布
}
generated quantities{
real xaste;
real log_lik;
xaste<- normal_rng(mu,sigma);           #予測分布
log_lik<-normal_log(x,mu,sigma);        #対数確率
}

実際、この.stanファイルを引数として受け取っているのは、G1mean()関数の中の次の行

fit<-stan(file=scr,data=dat,pars=par,seed=see,chains=cha,warmup=war,iter=ite)

で、これによって、推定結果がfitに返っています