INPUTしたらOUTPUT!

忘れっぽいんでメモっとく

第10回 「データ解析のための統計モデリング入門」 読書会に参加してきた

統計

ニコ生見直して復習していたので遅くなったけど以下メモ

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

(10/24 LTのタイトル間違えてたので修正。@piroyoung氏申し訳ございません）


第10章 「階層ベイズモデル」

前半


  • スライドの12枚目"対数尤度が小さくなったら、ある確率で対数尤度が大きくなるように更新する"とあるのがよく分からないな・・・
  • 階層ベイズモデル … 事前分布のパラメータにさらに事前分布が設定されている統計モデル
  • transformed parameters
    • Stanで他のパラメータから作るパラメータ
  • P.231の図10.4で中央値をプロットすると青丸のようになってしまった。。。
    • 平均値でプロットするとテキストのようになる
  • rは無事前情報分布でsは一様分布にするのはなぜか？
    • P.240に書いてある
    • sの事前分布としては逆ガンマ分布が使われることが多い
    • 逆ガンマ分布ではパラメータεに依存して事後分布が変わるので連続一様分布にした


後半


  • ベイズ統計モデルの設計で重要なこと
    • 事前分布の選択
  • よく使われる3種類の事前分布
    • 主観的な事前分布
      • 「私はこうなっていると思う!」という感じで事前分布をえいやと決めてしまう
        • 全然サイエンスではないので緑本では扱わない
      • 使わざるを得ない場合もある
        • 連続値の観測値に関する測定値の誤差が大きい場合など
        • 工学的なところでは使うこともある
    • 無情報事前分布
      • 事前分布が事後分布に影響を与えたくない場合に使用する
      • 統計モデリングではない
    • 階層事前分布
      • 一階層深くした事前分布は無情報事前分布にすることによって階層を深くしない
  • 事前分布の選び方
    • データ全体を大域的に説明する少数のパラメータ
      • 無情報事前分布
    • データの一部を説明する局所的パラメータ
      • 階層事前分布
  • 統計モデルと推定方法を区別する
    • 統計モデル : 階層ベイズモデル
    • 推定方法 : MCMCサンプリング
      • MCMCを使ったモデルというと混同していると言われる
  • 説明変数が連続量の場合、誤差を付けた方がよい？
    • 誤差項がないほうが良いと思う


LT

LaplacesDemonによるベイズ統計モデリング


  • 第43回R勉強会＠東京で聞いているのでメモ割愛


可視化で理解するマルコフ連鎖モンテカルロ法

超要約Stan Reference


  • increment_log_prob
    • 一番重要な関数
    • 対数の和を渡すことで事後分布をサンプルする
  • vectorize使うと速くなる
  • 自分にはビルの夜景に見えました
  • 7章と10章でモデルが異なるのでGLMMと階層ベイズモデルの比較が自分の中でできていない
    • 7章では{logit \left( q_i \right) = \beta_1 + \beta_2 x_i + r_i}
    • 10.3章では{logit \left( q_i \right) = \beta + r_i}
  • glmmMLで切片と個体差だけのモデルの作り方が分からなかったので階層ベイズモデルの方に葉数を入れてみた。
# Stanのコード
data {
  int<lower=0> N;
  int<lower=0> X[N];
  int<lower=0> Y[N];
}
parameters {
  real beta1;
  real beta2;
  real r[N];
  real<lower=0> sigma;
}
transformed parameters {
  real q[N];
  for(i in 1:N){
    q[i] <- inv_logit(beta1 + beta2*X[i] + r[i]); // 生存確率
  }
}
model {
  Y ~ binomial(8, q); // 二項分布
  beta1 ~ normal(0, 100);     // 無情報事前分布
  beta2 ~ normal(0, 100);     // 無情報事前分布
  r ~ normal(0, sigma);       // 階層事前分布
  sigma ~ uniform(0, 1.0e+4); // 無情報事前分布
}

# Rのコード
> library(rstan)
> d <- read.csv(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/glmm/data.csv"))
> d.dat <- list(N=dim(d)[1], X=d$x, Y=d$y)
> d.fit <- stan(file='kubo10.stan', data=d.dat, iter=10100, warmup=100, thin=10, chains=1)


なんかアラートが出たけどとりあえず動いた。

> d.fit
Inference for Stan model: kubo10.
1 chains, each with iter=10100; warmup=100; thin=10; 
post-warmup draws per chain=1000, total post-warmup draws=1000.

          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
beta1    -4.27    0.03 0.97   -6.30   -4.90   -4.19   -3.63   -2.44   951 1.00
beta2     1.02    0.01 0.23    0.61    0.86    1.01    1.17    1.51   801 1.00
r[1]     -2.01    0.06 1.79   -5.94   -3.06   -1.85   -0.77    0.94   885 1.00
r[2]     -0.08    0.04 1.19   -2.54   -0.86   -0.03    0.72    1.98   752 1.00
(以下略)


  • 中央値で見ると{\beta_1} は真の値-4に対して-4.19（GLMMでは-4.13）、{\beta_2} は真の値1に対して1.01（GLMMでは0.99）になっているのでそれなりに上手くいっているのかな？


  • 次回はいよいよ最終回！
  • 皆勤賞目指して頑張る