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

• ２章のコードを動かしてみる
  • G1mean()関数は、この教科書が提供しているstanを使う関数。以下のような設定で回るように作られている
    • 正規分布を仮定し
    • 正規分布のパラメタmu,sigmaをサンプリングしつつ
    • そのほかにxaste,log_likともサンプリングする
    • xasteは、そのときどきに選ばれたmu,sigmaのもとでの正規分布からの標本、log_likはそのときどきに選ばれたmu,sigmaのもとで、観測データxを観測する対数尤度

source('myfunc/myfunc.R')      #自作関数の読み込み
library(rstan)                 #パッケージrstanの呼び出し

#「牛丼データ」の入力
x<-c(76.5,83.9,87.9,70.8,84.6,85.1,79.6,79.8,79.7,78.0)

#正規分布に関する推測
#mu,sigmaの上下限を指定して実行
out <-G1mean(x,prior=T,mL=0, mH=1000, sL=0, sH=100)

• • 21000サンプリングして、初めの1000は捨てる設定で、それを５回繰り替えす条件で回す。その経過を表示

> out <-G1mean(x,prior=T,mL=0, mH=1000, sL=0, sH=100)

SAMPLING FOR MODEL 'G1meanPT' NOW (CHAIN 1).

Chain 1, Iteration:     1 / 21000 [  0%]  (Warmup)
Chain 1, Iteration:  1001 / 21000 [  4%]  (Sampling)
Chain 1, Iteration:  3100 / 21000 [ 14%]  (Sampling)
Chain 1, Iteration:  5200 / 21000 [ 24%]  (Sampling)
Chain 1, Iteration:  7300 / 21000 [ 34%]  (Sampling)
Chain 1, Iteration:  9400 / 21000 [ 44%]  (Sampling)
Chain 1, Iteration: 11500 / 21000 [ 54%]  (Sampling)
Chain 1, Iteration: 13600 / 21000 [ 64%]  (Sampling)
Chain 1, Iteration: 15700 / 21000 [ 74%]  (Sampling)
Chain 1, Iteration: 17800 / 21000 [ 84%]  (Sampling)
Chain 1, Iteration: 19900 / 21000 [ 94%]  (Sampling)
Chain 1, Iteration: 21000 / 21000 [100%]  (Sampling)# 
#  Elapsed Time: 0.031 seconds (Warm-up)
#                0.359 seconds (Sampling)
#                0.39 seconds (Total)
# 

SAMPLING FOR MODEL 'G1meanPT' NOW (CHAIN 2).

Chain 2, Iteration:     1 / 21000 [  0%]  (Warmup)
Chain 2, Iteration:  1001 / 21000 [  4%]  (Sampling)
Chain 2, Iteration:  3100 / 21000 [ 14%]  (Sampling)
Chain 2, Iteration:  5200 / 21000 [ 24%]  (Sampling)
Chain 2, Iteration:  7300 / 21000 [ 34%]  (Sampling)
Chain 2, Iteration:  9400 / 21000 [ 44%]  (Sampling)
Chain 2, Iteration: 11500 / 21000 [ 54%]  (Sampling)
Chain 2, Iteration: 13600 / 21000 [ 64%]  (Sampling)
Chain 2, Iteration: 15700 / 21000 [ 74%]  (Sampling)
Chain 2, Iteration: 17800 / 21000 [ 84%]  (Sampling)
Chain 2, Iteration: 19900 / 21000 [ 94%]  (Sampling)
Chain 2, Iteration: 21000 / 21000 [100%]  (Sampling)# 
#  Elapsed Time: 0.016 seconds (Warm-up)
#                0.405 seconds (Sampling)
#                0.421 seconds (Total)
# 

SAMPLING FOR MODEL 'G1meanPT' NOW (CHAIN 3).

Chain 3, Iteration:     1 / 21000 [  0%]  (Warmup)
Chain 3, Iteration:  1001 / 21000 [  4%]  (Sampling)
Chain 3, Iteration:  3100 / 21000 [ 14%]  (Sampling)
Chain 3, Iteration:  5200 / 21000 [ 24%]  (Sampling)
Chain 3, Iteration:  7300 / 21000 [ 34%]  (Sampling)
Chain 3, Iteration:  9400 / 21000 [ 44%]  (Sampling)
Chain 3, Iteration: 11500 / 21000 [ 54%]  (Sampling)
Chain 3, Iteration: 13600 / 21000 [ 64%]  (Sampling)
Chain 3, Iteration: 15700 / 21000 [ 74%]  (Sampling)
Chain 3, Iteration: 17800 / 21000 [ 84%]  (Sampling)
Chain 3, Iteration: 19900 / 21000 [ 94%]  (Sampling)
Chain 3, Iteration: 21000 / 21000 [100%]  (Sampling)# 
#  Elapsed Time: 0.016 seconds (Warm-up)
#                0.406 seconds (Sampling)
#                0.422 seconds (Total)
# 

SAMPLING FOR MODEL 'G1meanPT' NOW (CHAIN 4).

Chain 4, Iteration:     1 / 21000 [  0%]  (Warmup)
Chain 4, Iteration:  1001 / 21000 [  4%]  (Sampling)
Chain 4, Iteration:  3100 / 21000 [ 14%]  (Sampling)
Chain 4, Iteration:  5200 / 21000 [ 24%]  (Sampling)
Chain 4, Iteration:  7300 / 21000 [ 34%]  (Sampling)
Chain 4, Iteration:  9400 / 21000 [ 44%]  (Sampling)
Chain 4, Iteration: 11500 / 21000 [ 54%]  (Sampling)
Chain 4, Iteration: 13600 / 21000 [ 64%]  (Sampling)
Chain 4, Iteration: 15700 / 21000 [ 74%]  (Sampling)
Chain 4, Iteration: 17800 / 21000 [ 84%]  (Sampling)
Chain 4, Iteration: 19900 / 21000 [ 94%]  (Sampling)
Chain 4, Iteration: 21000 / 21000 [100%]  (Sampling)# 
#  Elapsed Time: 0.015 seconds (Warm-up)
#                0.406 seconds (Sampling)
#                0.421 seconds (Total)
# 

SAMPLING FOR MODEL 'G1meanPT' NOW (CHAIN 5).

Chain 5, Iteration:     1 / 21000 [  0%]  (Warmup)
Chain 5, Iteration:  1001 / 21000 [  4%]  (Sampling)
Chain 5, Iteration:  3100 / 21000 [ 14%]  (Sampling)
Chain 5, Iteration:  5200 / 21000 [ 24%]  (Sampling)
Chain 5, Iteration:  7300 / 21000 [ 34%]  (Sampling)
Chain 5, Iteration:  9400 / 21000 [ 44%]  (Sampling)
Chain 5, Iteration: 11500 / 21000 [ 54%]  (Sampling)
Chain 5, Iteration: 13600 / 21000 [ 64%]  (Sampling)
Chain 5, Iteration: 15700 / 21000 [ 74%]  (Sampling)
Chain 5, Iteration: 17800 / 21000 [ 84%]  (Sampling)
Chain 5, Iteration: 19900 / 21000 [ 94%]  (Sampling)
Chain 5, Iteration: 21000 / 21000 [100%]  (Sampling)# 
#  Elapsed Time: 0.031 seconds (Warm-up)
#                0.39 seconds (Sampling)
#                0.421 seconds (Total)
# 

• • outオブジェクトにすべてが格納してある
  • G1mean()関数の中身をみると

out<-list(fit=fit,par=par,prob=prob,mu=mu,sigma=sigma,xaste=xaste,log_lik=log_lik)
return(invisible(out))

• • とあるので、out$fit (out1)はstanの出力そのもの、mu(out4),sigma(out5),xaste(out6,log_lik(out7)は４種類の記録対象の事後標本であることがわかる

str(out)

• • mu,sigma,xaste,log_likをプロットすると100,000標本があることがわかるので、初期1000を捨てた後の(21000-1000)*5の事後分布標本であることがわかる

par(mfcol=c(2,2))
plot(out[[4]])
plot(out[[5]])
plot(out[[6]])
plot(out[[7]])
par(mfcol=c(1,1))

• • print()関数が用意されているので、教科書が指定した「知りたいこと」はプリントアウトできる
  • では、stan()関数が残して来れば、事後分布サンプルから同じことをすることはできるだろうか

summary(out[[4]])
quantile(out[[4]],c(0.025,0.05,0.5,0.95,0.975))
summary(out[[5]])
quantile(out[[5]],c(0.025,0.05,0.5,0.95,0.975))
summary(out[[5]]^2)
summary(out[[5]]/out[[4]])

• • 観測値が80未満になる確率

# 事後分布標本から計算
length(which(out[[6]]<80))/length(out[[6]])
# 正規分布としての事後分布標本から、条件を満たすサンプリングが起きる確率の平均
q <- 0
for(i in 1:length(out[[4]])){
   q <- q + pnorm(80,out[[4]][i],out[[5]][i])/length(out[[4]])
}
# または
q2 <- mean(pnorm(80,out[[4]],out[[5]]))

• rstanをもっと自由に動かしたいなら→こちら