準備:ぱらぱらめくる『統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―』
はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―
- 作者: 豊田秀樹
- 出版社/メーカー: 朝倉書店
- 発売日: 2016/06/02
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (11件) を見る
- 準備
- こちらに書籍用の諸ファイルがあります
- ダウンロードしたzipファイルを展開すると、
- 3つのフォルダ
- myfunc
- srcA
- stan
- 2つのファイル
- Readme.txt
- RstanStan.pdf
- 3つのフォルダ
- があります
- Readme.txtは、このフォルダとその中の諸ファイルを使うための準備について書いてあります
- 必要なものは
- R本体
- rstanというパッケージ
- 必要なものは
install.packages("rstan")
- その準備が揃ったところで
- Working directoryを、このzipを展開したところ(statという名前のディレクトリ)に設定します
- そのstatディレクトリを指定します。Desktopに展開したなら、たとえば
setwd("C:/Users/ryamada/Desktop/stat")
-
- のように。
- Rstanのインストールについてはこんなページもあります
実行:ぱらぱらめくる『統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―』
- 実行は、srcAディレクトリのファイルをコピーペーストして実行するだけ
- 6章分
source('myfunc/myfunc.R')
- 冒頭でmyfuncディレクトリのmyfunc.Rソースを読み込みます
- このmyfunc.Rは、myfuncディレクトリ内のそのほかのいろいろなxxx.Rソースを読み込みます
- 第1章用のファイル"part1.R"はstanを使うより前
- 第2章用から、stanを使います。それは第2章用ファイル"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に返っています