ベイズ統計の学習記録（演習#6-3-3） - Fake Data-Scientistの学習日記

章末問題6.3-3)
題意より、昨年大会のボーダーラインを0.90分位数のEAP推定値を以下で求めます。

$\begin{align} \displaystyle \xi_{0.90}^{(t)} = g(\mu^{(t)}, \sigma^{(t)}) = \mu^{(t)} + z_{0.90}\sigma^{(t)} = \mu^{(t)} + 1.282\sigma^{(t)} \end{align}$

続いて、B子さんの実績(平均87点、標準偏差5点)より、B子さんが本選に進める確率を定式化すると、

$\displaystyle \begin{align} p_{1-\Psi(\xi_{0.90})}^{(t)} = g(\xi_{0.90}^{(t)}) = 1 - \Psi(\xi_{0.90}^{(t)}) \end{align}$

のようになります。

ここで、normal_cdf関数を用いれば、 $1-\Psi(z)$ に標準化( $z = \frac{\xi_{0.90}^{(t)}-87}{5}$ )された正規累積分布関数の確率を返してくれます。

実際に使用したStanコードとRコードは以下の通りです。

Stanコード

#ex_6-3-3.stan

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

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

#パラメータの数値変換
#0.90分位数を平均と標準偏差より算出
transformed parameters {
  real xi;
  xi = mu + 1.282*sigma;
}

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

#生成量の計算
#0.90分位数を超える確率を計算
generated quantities{
  real<lower=0, upper=1> prob_over;
  prob_over = 1- normal_cdf(xi , 87, 5);
}

Rコード

library(rstan)

#データセットの作成
karaoke_score <- c(75, 82, 77, 86, 98, 91, 85, 84, 82, 79,
                88, 79, 84, 87, 69, 93, 84, 82, 89, 78,
                90, 74, 75, 84, 89, 81, 74, 79, 82, 84)

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

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

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

#Stanのmodelファイルをコンパイル
fit_stan <- sampling(
  stanmodel,
  data = data,
  par=pars,
  seed = 1234,
  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-3.
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
xi        91.297   0.019 1.714 88.348 90.100 91.155 92.337 95.024  7962    1
prob_over  0.208   0.001 0.089  0.054  0.143  0.203  0.268  0.394  8772    1