２−５．rstanを自由に動かす　ぱらぱらめくる『はじめての統計データ分析』

教科書のコードをコピーペーストすれば、MCMCによる事後分布推定はできるけれど、条件を変えにくい
少し勉強して自由に動かしてみる
rstanパッケージのstan()関数に必要なものを渡せばよい
- どうしても渡さないといけない引数はfile、dataの二つ
- fileはどういう分布をモデルにしているか、パラメタの名前、記録するべき値、などMCMC実行条件を記載したファイル(へのパス)
- dataは観測データを含むlist。値ベクトル(値行列)を要素とするlistで、要素数も明示するのがルールらしい
- やってみる

# 教科書のG1mean()関数が用いているモデル指定ファイル
file. <- "stan/G1meanPF.stan"
# 観測データである値ベクトルxをlistにする
x.list <- list(n=length(x),x=x)
fit2 <- stan(file=file.,data=x.list)
# 事後分布標本はextract()関数を用いて「取り出す」
ext.samples <- extract(fit2)
plot(ext.samples$mu)
par(mfcol=c(2,2))
plot(ext.samples[[1]])
plot(ext.samples[[2]])
plot(ext.samples[[3]])
plot(ext.samples[[4]])
par(mfcol=c(1,1))

後は"hoge.stan"ファイルの書き方がわかればよい。"G1meanPF.stan"の中身を見てみる
data,paarmeters,transformed parameters,model,generated quatitiesという要素からなることがわかる
dataは実観測データ
prametersはmodelのパラメタ。事後分布推定の対象
transfromed parametersは(おそらく、parametersから一意に決まるパラメタ
modelは、観測データが、parameterからどのように確率的に決まるかを規定する。parameterがどのような分布に従うか(事前分布を持つか)も指定する
generated quantitiesはいわゆる生成量。事後分布サンプリングされるべきもの
概念がわかったら、あとは、練習するのみ
- modelにて、観測データの尤度に相当するものを計算する必要がある。~を使って理論分布から計算することができる場合もあれば、そのような簡便記法がない場合もある
こちらやこちらなどが参考になる

data {
int<lower=0> n;                                  #データ数
real<lower=0> x[n];                              #データ
}
parameters {
real                    mu;             #平均(範囲指定)
real<lower=0>           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);        #対数確率
}