読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRでベイズ統計モデリングします. たまに書評.

「はじめての統計データ分析」 豊田秀樹のメモ

書評 Stan

あとがきと6章のあとにあるQ&Aの節が熱い思いに満ちていてオススメです。2.7節「論文・レポートでの報告文例」もユニークです。学生思いの教育者としての一面を垣間見た気がします。

あとがきに書いてあるように、たしかに初級向けの授業で伝統的な統計学と検定のラッシュを学び、中級以上向けの授業でベイズ統計モデリングを学ぶとしたら、内容の一貫性が乏しく、学ぶ側は(教える側も)違和感を覚えるかもしれません。その点、この本ではt検定に相当するような簡単なものから一貫してベイズ統計です。

また、各例題に対してリサーチクエスチョン(RQ)をきちんと設け、それに対してMCMCサンプルを使った生成量と予測分布を用いてシンプルに回答していくスタイルは、分かりやすくて読みやすいです。ベイズ統計の長所と思います。

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

しかし、予測分布に関しては一部分かりにくい箇所があります。具体的に言うと、2章以降の各章において確率変数である新しく取得するデータ  x^{*} が従う分布として、「事後予測分布」(2.3.1節)と「条件付き予測分布」(2.3.2節)の二つが大きな区別なく混在しているところです。

前者の事後予測分布はベイズの枠組みでは以下になります。

 p(x^{*} \mid X)=\int f(x^{*}\mid \theta) p(\theta \mid X) d \theta

ここで、  Xはデータ、  \thetaはパラメータ、  f(x^{*}\mid \theta)はモデルの分布、 p(\theta \mid X)は事後分布です。後者の「条件付き予測分布」は以下です。

 f(x^{*}\mid  \hat{\theta})

この \hat{\theta}には「最尤推定値」・「MAP推定値」・「事後平均の値(EAP)」といった分布の代表値を代入して使うことが多いです。渡辺洋氏の「ベイズ統計学入門」のp.69にもちらっと「条件付き予測分布」という単語は出てきますが、これは \hat{\theta}としてその3つの値を代入すること想定した書き方と思います。渡辺澄夫氏の「ベイズ統計の理論と方法」のp.16では、これらの推定値を代入して予測分布を構成する方法を、それぞれ「最尤推測」・「事後確率最大化推測」・「平均プラグイン推測」と呼んでいます。その本の後の章では、事後分布が(多変量)正規分布で近似できない場合には、これらで構成した予測分布は事後予測分布よりも予測の性能が汎化損失の観点からよくないことが示されています。なお、汎化損失は情報量規準のWAICと密接に関わっています。また「はじめての統計データ分析」に出てくるような簡単なモデルでは、事後分布は(多変量)正規分布で十分によく近似できるので、どの予測分布を使っても結果はほとんど変わらないことも書かかれています。 statmodeling.hatenablog.com しかし、「はじめての統計データ分析」ではそのいずれの構成方法とも異なり、MCMCサンプルの各値を代入して“予測分布”をMCMCサンプルの数だけ作って、その各々の“予測分布”の分位点(例えば25%点)のMCMCサンプルを作成し、これを“条件付き予測分布からの25%点の事後分布(分布の分布ですからメタ分布です)”(p.51)と呼んで、2章から6章で頻繁に使っています。この“メタ分布”は一般的に使われていないので注意が必要です。理由を以下に3つ挙げます。

“予測分布”が分布関数とは言えない

“条件付き予測分布からのX%点の事後分布”と呼んでいますが、新たに取得するデータ x^{*}の従う予測分布という分布関数がただ一つ存在するわけではありません。これを果たして“ x^{*}の予測分布”と呼んでいいのかよく分かりません。使用例が掲載されている論文もないとのことで、個人的には理論的な妥当性があるようには思えませんでした。もしこの“メタ分布”を使用した解析を統計関連の学会等で発表する場合には、発表の中できちんと説明する必要があるでしょう。説明がないと全く意思疎通できずに無視されるか攻撃される可能性があります。

予測の性能が悪い

“条件付き予測分布からのX%点の事後分布”はMCMCサンプルの各値を代入した f(x^{*}\mid \hat{\theta})がもとになっているため、事後分布が(多変量)正規分布でよく近似できない場合には事後予測分布との乖離が大きく、汎化損失の観点から予測の性能が悪くなります。

例えば、モデル式の一部が x ~ Normal(mu, 2) である場合を考えましょう。そして、muの事後分布が図中のピンクの線で得られたとしましょう。この場合にxの事後予測分布は緑の線になります。その95%点は図中の実線です。一方、“条件付き予測分布からの95%点の事後分布”はピンクの線を平行移動した青の線になります。その青の分布のMAP推定値と事後平均(図中の点線)は事後予測分布の95%点とは大きく離れてしまっています。

fig.png

「はじめての統計データ分析」では簡単なモデルに終始しているため、事後分布が(多変量)正規分布で十分によく近似できます。そのため、“条件付き予測分布からのX%点の事後分布”のMAP推定値や事後平均は事後予測分布のX%点とほぼ一致し、区間が求められる分だけ、“条件付き予測分布からのX%点の事後分布”の方が得に思えました。しかし、さらにベイズ統計モデリングをすすめれば事後分布が(多変量)正規分布で近似できない場合が普通です。一般的には理論的な妥当性のない上に予測の性能が悪い区間を使うより、素直に事後予測分布を使った方がよいでしょう。なお、データの多少を反映して区間の幅が変わることを見たければ、パラメータの事後分布や事後予測分布の区間の幅が変わることを見ればよいと思います。

正規分布以外の場合に算出が面倒

「はじめての統計データ分析」では正規分布を用いる例題多く、正規分布のX%点は平均パラメータと標準偏差パラメータから簡単に算出できるため、“条件付き予測分布からのX%点の事後分布”は求めやすいです。しかし、これがガンマ分布のX%点だったらどうでしょうか。累積分布関数の逆関数を使う必要があるため、Stan内では簡単にできず、R側で各MCMCサンプルに対してqgamma関数を使って求める必要があります。ベイズ統計モデリングをすすめれば正規分布以外を使うことも多いので、これは面倒です。

まとめ

「はじめての統計データ分析」は“メタ分布”の箇所以外はとても良い本だと思います。“条件付き予測分布からのX%点の事後分布”については上記3点のデメリットをよく理解して相応の覚悟があれば使ってもよいと思いますが、現実のデータ分析でベイズ統計モデリングしたい人には使わないことを強く推奨したいです。

ちなみに、事後予測分布は「MCMCサンプルを持っておく必要がある」「計算が大変」と言う人もいます。前者はたしかにその通りです。しかし、“条件付き予測分布からの25%点の事後分布”もMCMCサンプルを持っておく必要があります。後者に関してはStanではgenerated quantitiesブロック内で分布名_rng関数を使うことで、事後予測分布からのMCMCサンプルが簡単に得られます(Stan 2.9.0のマニュアルの37.3. Pseudorandom Number Generatorsの節を参照)。このブロックの計算は推定計算の100倍~1000倍ぐらい早い様子なので全く問題になりません。