ベイズ統計の学習記録（演習#6-3-1）

ベイズ統計の学習にあたり、以下の書籍で解いたプログラミング演習を備忘録の意味も込め、掲載しています。*1

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

章末問題6.3-1)
正規分布を用いて、12週の売上平均を推測する問題となっています。
実際に使用したStanコードとRコードは以下の通りです。

Stanコード

#ex_6-3-1.stan

#データサイズと確率変数Yを定義
data{
  int N;
  real Y[N];
}

#正規分布の平均と分散を設定
parameters{
  real mu;
  real<lower=0> sigma;
}

#確率変数Yを平均mu, 標準偏差sigmaの正規分布に従う確率変数Yを生成
model{
  for(i in 1:N){
    Y[i] ~ normal(mu, sigma);
  }
}

#生成量の計算
#delta_over: 12週の平均売上が7万円を超える確率を計算
generated quantities{
  real delta_over;
  real delta_over2;
  delta_over = step(mu - 70000);
  delta_over2 = step(mu - 75000);
}

Rコード

library(rstan)

#データセットの作成
sales <- c(76230, 73550, 80750, 71500, 75420, 74840,
           71580, 76920, 68450, 76990, 64070, 76200)

#Stanで読み込むデータに整形
data <- list(N=length(sales), Y=sales)

#Stanのmodelファイルをコンパイル
stanmodel <- stan_model(file='./ex_6-3-1.stan')

#表示するパラメータの設定
pars <- c("mu", "delta_over", "delta_over2")

#Stanのmodelファイルをコンパイル
fit_stan <- sampling(
  stanmodel,
  data = data,
  par=pars,
  init=function(){
    list(mu=runif(1, 0,1000), sigma=runif(1, 0, 5))
  },
  seed = 1234,
  #11000回のサンプリング、ワームアップは1000回
  chains=1, iter=11000, warmup=1000
)

#結果の表示
#有効桁数は3と設定
dig <- 3
print(fit_stan,pars=pars, digits_summary=dig)

#traceplotでMCMCサンプルが収束しているどうかを確認
traceplot(fit_stan, pars=pars, alpha=0.5)

結果はこんな感じでした。

Inference for Stan model: ex_6-3-1.
1 chains, each with iter=11000; warmup=1000; thin=1; 
post-warmup draws per chain=10000, total post-warmup draws=10000.

                 mean se_mean       sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
mu          73873.476  21.176 1517.223 70825.24 72927.85 73877.61 74795.38 76959.06  5133    1
delta_over      0.992   0.001    0.091     1.00     1.00     1.00     1.00     1.00  4970    1
delta_over2     0.207   0.005    0.405     0.00     0.00     0.00     0.00     1.00  5862    1

*1:あくまでの個人の見解ですので、間違いのある可能性もございます。その場合はご了承ください。