ベイズモデリングの実践編_#1-2
「基礎からのベイズ統計学」と並行して、「ベイズモデリングの実践編」も学習中です。
徐々に、Stanの使い方やReferenceの見方をキャッチアップできてきました。
- 作者: 豊田秀樹
- 出版社/メーカー: 朝倉書店
- 発売日: 2017/01/25
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (2件) を見る
早速ですが、実装した教科書の例題(分析例: 1-2)を載せていきます。
ガンベル分布を利用して、極限統計を解く問題です。
使用したStanコードとRコードは以下となります。
Stanコード
#ana1-2.stan data{ int<lower=0> N; real<lower=0> Y[N]; } parameters{ #RQ.1 最頻値と確信区間 real mu; real<lower=0> sigma; } model{ for(i in 1:N){ Y[i] ~ gumbel(mu, sigma); } } generated quantities{ real<lower=0> s; real p; real x; real<lower=0, upper=1> U; real inv_r; real y; real<lower=0, upper=1> Uy; #RQ.2 各都市の最長記録の散らばり(標準偏差) s = pi() * sigma * inv_sqrt(6); #RQ.3 実現値の確率 p = exp(-exp((mu-8.95)/sigma)); #Rq.4 ガンベル分布の99%点の再現レベル x = mu - sigma * log(-log(0.99)); #Rq.4 世界記録の再現期間が100年より長い確率 U = x < 8.95 ? 1 : 0; #RQ.5 再現期間の必要年数 inv_r = inv(1 - p); #RQ.6 推定パタメータを用いたガンベル分布の予測値 y = gumbel_rng(mu, sigma); #翌年に世界記録が更新される確率 Uy = y > 8.95 ? 1 : 0; }
Rコード
library(rstan) #各年の最長記録データ max_record <- c(8.95, 8.68, 8.70, 8.74, 8.71, 8.58, 8.63, 8.48, 8.60, 8.65, 8.41, 8.52, 8.53, 8.60, 8.60, 8.56, 8.66, 8.73, 8.74, 8.47, 8.54, 8.35, 8.56, 8.51, 8.52) #読込データの設定 data <- list(N=length(max_record), Y=max_record) #モデルのコンパイル stanmodel <- stan_model(file='./ana1-2.stan') #表示パラメータの設定 pars <- c("mu","sigma", "s", "p", "x", "U", "inv_r", "y", "Uy") #サンプリング fit_stan <- sampling( stanmodel, data = data, par=pars, seed = 1234, chains=4, iter=20000, warmup=5000 ) #結果確認 dig <- 3 print(fit_stan,pars=pars, digits_summary=dig) #トレースプロット traceplot(fit_stan, pars=pars, alpha=0.5)
結果はこんな感じでした。
Inference for Stan model: ana1-2. 4 chains, each with iter=20000; warmup=5000; thin=1; post-warmup draws per chain=15000, total post-warmup draws=60000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat mu 8.542 0.000 0.026 8.492 8.525 8.542 8.559 8.594 30623 1 sigma 0.121 0.000 0.019 0.090 0.107 0.119 0.133 0.166 28628 1 s 0.156 0.000 0.025 0.115 0.138 0.153 0.170 0.213 28628 1 p 0.963 0.000 0.023 0.905 0.952 0.968 0.979 0.991 24923 1 x 9.100 0.001 0.100 8.935 9.030 9.090 9.159 9.326 26889 1 U 0.041 0.001 0.198 0.000 0.000 0.000 0.000 1.000 36026 1 inv_r 39.402 0.171 30.153 10.497 20.848 31.157 48.012 117.175 31026 1 y 8.613 0.001 0.161 8.378 8.502 8.586 8.694 9.007 55916 1 Uy 0.038 0.001 0.191 0.000 0.000 0.000 0.000 1.000 60000 1
最後にとのEAP推定値を用いて、ガンベル分布の密度関数をプロットします。
library(evd) x <- 8:10 dgumbel(x, loc=8.542, scale=0.121) curve(dgumbel(x, loc=8.542, scale=0.121), 8, 10, type="l", ylab="f(x)")