Fake Data-Scientistの学習日記

似非データサイエンティストの学習記録です。

ベイズモデリングの実践編_#1-2

「基礎からのベイズ統計学」と並行して、「ベイズモデリングの実践編」も学習中です。
徐々に、Stanの使い方やReferenceの見方をキャッチアップできてきました。

実践 ベイズモデリング -解析技法と認知モデル-

実践 ベイズモデリング -解析技法と認知モデル-


早速ですが、実装した教科書の例題(分析例: 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


最後に \mu \sigmaEAP推定値を用いて、ガンベル分布の密度関数をプロットします。

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)")

f:id:FAKED-SCIENTIST:20180504210928p:plain