下校時刻

メモ、読んだ本のまとめメモ、考えたことのメモ、その他のメモなどが書いてあるブログです(twitter: @hoture6)

「StanとRでベイズ統計モデリング」まとめ

https://www.amazon.co.jp/Stan%E3%81%A8R%E3%81%A7%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E3%83%A2%E3%83%87%E3%83%AA%E3%83%B3%E3%82%B0-Wonderful-R-%E6%9D%BE%E6%B5%A6-%E5%81%A5%E5%A4%AA%E9%83%8E/dp/4320112423/ref=pd_sim_14_8?_encoding=UTF8&psc=1&refRID=3GZ40XY6NJA51QHPBAGYwww.amazon.co.jp

まとめます。

読んでいる途中から、自分が文法的に耐えられるのはRよりPythonだしStanよりPyMCだということに気づいて若干やる気が失せたので、後半は適当に読んでしまいました。しかし、統計モデリングをする上での作法・考え方のようなものがちゃんとした事例をもとにわかりやすく述べられていて、かなり良い本だと思います。また、統計の本によくある確率分布を羅列する章がこの本にもあるのですが、いままで読んだ確率分布を羅列する章の中でこの本の確率分布を羅列する章がもっともやる気が失せずに読めてかつ役立つ感じでした(具体的な使い所が書いてあるのと性質が表にまとまっているところが良かった)。今後Stanを使った方が良い場面もあると思うので、Stanの文法や基本的な使い方がわかったもの良かったです。


Ch. 1 「統計モデリングとStanの概要」

統計モデリングとは

  • モデル:不必要な性質を無視して本質だけを取り上げたもの
  • 数理モデル:数式を使ったモデル
  • 確率モデル:確率分布・パラメータとその関係式からなる数理モデル

統計モデリングとは、確率モデルをデータに当てはめることで現象の理解と予測を促す営み。データと確率分布を入力とし、パラメータを推測する。

機械学習・古典的な統計分析との違い

  • 機械学習:数式に知識を組み込む作業がないため、人間にとってよくわからん現象でもパフォーマンスが高いが、結果の解釈が難しい
  • 古典的な統計分析:結果の解釈はしやすいが、予測に役立てるのは難しい
  • 統計モデリング:結果の解釈もしやすく、予測性能も高い

Stanの特徴

Hamiltonianモンテカルロの一つの実装であるNUTSを使っているため、ステップ間の相関が弱い。また変分ベイズ法の一つの実装であるADVIも使える(近似だけど速い)。

Ch. 2 「ベイズ推定の復習」

最尤推定ベイズ推定

最尤推定では、パラメータ$\theta$は固定値だと考えて、尤度$L(\theta) = p(y = Y|\theta)$を最大にする$\theta$を求める。ここで、$L(\theta)$は確率ではないので1よりも大きくなり得ることに注意(離散分布だと確率になる?)。

ベイズ推定では、パラメータ$\theta$は確率変数だと考えて、事後分布$p(\theta|Y) \propto p(Y|\theta)p(\theta)$を求める。無情報事前分布を用いた場合のMAP推定値(事後分布が最大になる点)は最尤推定値と一致する。

MCMC

MCMCは事後分布を求めるための手法である。重要な用語として以下がある。

  • chain:初期値と乱数の種を固定したサンプル列
  • warm up / burn in:初期値の影響を受けているためサンプル列から捨てるステップ数
  • thinning:サンプルの相関を下げるために何個かおきにサンプルすること
  • convergence:ステップ数に対して事後分布が変化しなくなり、かつすべてのchainで同じ分布になっていること

Ch. 3 「統計モデリングをはじめる前に」

統計モデリングの手順

  • 解析の目的を設定
  • データの分布の確認
  • メカニズムの想像
  • モデル式の記述
  • シミュレーション
  • MCMC
  • グラフによるモデルのチェック

背景知識の役割

統計モデリングでは、(一部の)インプットとアウトプットが与えられて、その間のメカニズムがわかっていない現象を解析することが多い。メカニズムを表すのに何かしらモデルを設定する必要があるが、現実の問題では「真のモデル」は知り得ないため、「解釈・納得しやすく、予測がちゃんとできそうなモデル」を背景知識をもとに設定する。

Ch. 4 「StanとRStanをはじめよう」

Stanのコード

data {
    入力データX、出力データY、その他のデータを列挙
}

parameters {
    modelに使うパラメータ(=推定したいパラメータ)を列挙
}

transformed parameters {
    parametersを便利に変換したもの
}

model {
    Yはどの分布に従うか、そのパラメータは何か
}

generated quantities {
    サンプルしたい確率変数
}

Ch. 5 「基本的な回帰とモデルのチェック」

重回帰

説明変数が複数の回帰。説明変数はスケーリングして同程度の大きさに変換した方がよい。また、結果の確認の際には「横軸:実測値、縦軸:予測値」のグラフをベイズ予測区間とともに書くと良い(!)。

ロジスティック回帰

ポアソン回帰

Ch. 7 「回帰分析の悩みどころ」

変数変換

結果の解釈を難しくする変数変換はしない方がよいが、対数変換は有用なことが多い。回帰係数の意味としては、

  • 普通に回帰:Xを+xしたとき時Yは何増えるか
  • 対数をとって回帰:Xをa倍した時Yは何倍になるか

となる。もちろん、どちらのモデルが正しいとか間違っているということはない。

交互作用

重回帰分析において説明変数同士の掛け算の項を考えること。解釈が難しくなるので、データをプロットしてみて明らかに交互作用が確認できる場合以外は考慮しない方がいい。

多重共線性

重回帰分析において説明変数間の相関が高すぎて回帰係数が収束しないこと。例えば、$\mu[n] = a A[n] + b B[n]$において、$A[n] \simeq B[n]$だとすると、近似的に$\mu[n] = (a + b) A[n]$となってしまうので、$a, b$は一意に決まらなくなる。これを防ぐために、相関が明らかな場合は片方の変数を捨てる。

交絡

モデルの外側に、応答変数と説明変数の両方に影響を与える変数があること。パス解析をすることで確認できる。

外れ値への対処

  • 外れ値を取り除いて解析をする
  • 裾の長い(まれに大きな値を生成する)確率分布を使って解析をする

Ch. 8 「階層モデル」

階層モデルの導入

ある業種のいくつかの似たような会社について、年齢と収入の関係を予測する。

  1. 会社間の差を無視して回帰 → 会社ごとの傾向の違いを考慮できない
  2. 会社ごとに回帰 → データ数が少なくなるので外れ値に弱いと同時に、年齢と収入の関係について一般的な知見が得られない。
  3. 階層モデルで回帰 → いい。

複数階層

さらに他の業種のいくつかの会社も同時にモデリングしたい場合、もう一階層増える。

Ch. 10 「収束しない場合の対処法」

パラメータの識別可能性

同一の事後分布を与えるパラメータの値が複数あるとき、そのパラメータは識別不可能であるという。以下のような場合に生じる。

  • データが各個人に一つしかないのに個人差を考慮してしまう場合(ノイズと個人差の区別がつかない)
  • ラベルスイッチング
  • softmax関数への入力の全体的なシフト → 一つの成分を固定する!

弱情報事前分布

少ないデータを複雑にモデリングすると、無情報事前分布では収束しない場合がある。弱情報事前分布を使うとよい(ただし、最初から使うのではなく無情報でうまくいかないとき使ってみるという感じ)。

  • [a, b]の範囲の可能性が高いパラメータ:$\mu = (a + b) / 2$、$\sigma = (b - a) / 2$の正規分布
  • 回帰係数:大きすぎる値を避けるためt分布
  • 正の値をとるパラメータ:半t分布や指数分布
  • 確率:ベータ分布やディリクレ分布