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

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

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

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

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

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

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

ここで、 はデータ、 はパラメータ、 はモデルの分布、は事後分布です。後者の「条件付き予測分布」は以下です。

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

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

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

予測の性能が悪い

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

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

「はじめての統計データ分析」では簡単なモデルに終始しているため、事後分布が（多変量）正規分布で十分によく近似できます。そのため、“条件付き予測分布からの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倍ぐらい早い様子なので全く問題になりません。

Michael Betancourt's Stan Lectureを開催しました

ドワンゴさんに会場提供していただき、2016/6/4にMichael Betancourt's Stan Lectureを開催しました。実はStanの勉強会というのはこれがはじめてではなく、約2年ほど前に催されていたBUGS/Stan勉強会がもととなっています。

また今回はニコ生で放送したのですが、400人以上の視聴者がいたのは驚きました。この内容で、しかも英語なのに。欲を言うともっと会場にも来てほしかったです。当日の進行は少し前倒しになり、そのために閲覧できなかった方には申し訳ありませんでした。動画を公開しますのでお許しください。

以下では内容を簡単に紹介します。

Hiroki ITO "Dealing with latent discrete parameters in Stan"

一人目は頼み込んで北海道から（自腹で）来てくれました伊東宏樹さんです。動画は次の僕の分と合わせて以下になります。

スライドは以下です。

Dealing with latent discrete parameters in Stan from Hiroki Itô

講演のもととなっているのは、伊東さんが一部翻訳された以下の本です。

BUGSで学ぶ階層モデリング入門: 個体群のベイズ解析

講演の内容は前半が「data augmentation」、後半が「隠れマルコフモデル」です。ここでは「data augmentation」の説明を簡単にします。

ある調査地で生物の個体をカウントすることを考えます。1回目は30匹いたとしましょう。その30匹に印をつけてまた放します。別の日に同じ調査地でまた生物の個体をカウントします。前に印をつけた個体も含まれていることでしょう。このようにカウントを複数回行なったデータがあるとします。その調査地に存在しないかつ複数回の観測で一回も発見されなかった仮想的な個体のデータを十分に加え（data augmentation）、「その調査地にそもそも個体が存在する確率」と「観測時に（カウント時に）発見する確率」を切り分けてモデリングします。すると、それらの確率を推定できるだけでなく、その調査地における総個体数も見事に推定できます。この考え方はとても柔軟で応用範囲が広く、上記の本の見どころの一つとなっています。個人的には上記の本で一番感動した部分です。

上記の本はBUGSの本ですが、伊東さんがすべてStanコードに翻訳し、Stanの公式リポジトリにpushされています。 github.com 上記リンク先のCh.06のフォルダが「data augmentation」の初出の章になっています。

Kentaro Matsuura "Replica exchange MCMC with Stan and R"

スライドは以下です。

Replica exchange MCMC from . .

詳細のブログ記事は以下です。 statmodeling.hatenablog.com 質問の答えをここで補足します。simulated annealingでは最適値を求めますが、Replica exchange MCMCでは分布を求めることができます。それが大きな違いです。

Michael Betancourt "Scalable Bayesian Inference with Hamiltonian Monte Carlo"

動画は以下になります。

前半のTall data, Wide dataの概論から多数のパラメータをもつモデルでも効率的にサンプリングできるHMCについて話をすすめていきます。HMCでは仮想的に運動量を導入するのですが、この部分は参加者のxiangzeさんの記事が詳しいです。 xiangze.hatenablog.com NUTSはHMCの説明で出てきた各パラメータを自動で調節するアルゴリズムです。その詳細はHMCで脱落するぐらい人達にはムズすぎるのでカットされた模様です。また、最後の方にReplica exchange MCMCにあわせてかAdiabatic Monte Carlo (arXiv)を紹介していただきました。

Michael Betancourt "Some Bayesian Modeling Techniques in Stan"

動画は以下になります。

線形モデル・GLMから階層モデルまで丁寧な説明でした。後半では「centered vs. non-centered parameterization」の話が面白かったです。centeredは普通の階層モデルで、non-centeredはスケールを分離するような再パラメータ化をしたモデルです。詳しくは動画やStanのJSS論文を見てください。

StanのJSS論文から一部引用します。

Betancourt and Girolami (2013) provide a detailed analysis of the benefits of centered vs. non-centered parameterizations of hierarchical models in Euclidean and Riemannian HMC, with the conclusion that non-centered are best when data are sparse and centered when the data strongly identifies coefficients.

しかしながらデータのスパースネスを判定する方法は一般には難しいので、推奨は以下のようになるでしょう（JSS論文からの引用）。

We recommend starting with the more natural centered parameterization and moving to the non-centered parameterization if sampling of the coefficients does not mix well.

また、とりあえず全部centeredにしておいて、ある階層だけ・あるグループだけ部分的にnon-centeredにするとよい場合も考えられるそうです。とても勉強になりました。

最後に、会場＆ニコ生を手配していただいたSAMさんとドワンゴさん、司会や受付をしてくれた若者たち、伊東さん、Michael Betancourtさん、本当にありがとうございました！

Report of Michael Betancourt's Stan Lecture

We are really happy to hold Michael Betancourt's Stan Lecture on June 4 at DWANGO. To tell the truth, this is not the first Stan meeting in Japan. Three BUGS/Stan meetings were held about 2 years ago.

DWANGO is a very famous company for its video service (nicovideo), and the following lectures were broadcast using it. I'm surprised that there were 400+ listeners even contents were English and specialized.

I'd like to introduce the contents.

Hiroki ITO "Dealing with latent discrete parameters in Stan"

That movie includes the next presentation. The slides are here.

Dealing with latent discrete parameters in Stan from Hiroki Itô

This lecture was based on the book Hiroki translated:

BUGSで学ぶ階層モデリング入門: 個体群のベイズ解析

The first half of his presentation was on "data augmentation", and the last half was on "Hidden Markov Model". Data augmentation is an amazing technic to estimate the probability of existence and the probability of observation separately by adding virtual individuals' data. Furthermore, the number of existing individuals can be also estimated. This technic is very flexible and widely applicable. I think data augmentation is one of the best parts in the above book.

The above book is about BUGS, but Hiroki has translated all of them into Stan, and pushed them to the official repository of Stan. github.com Ch.06 is the first directory that contains source codes of data augmentation.

Kentaro Matsuura "Replica exchange MCMC with Stan and R"

The slides are here.

Replica exchange MCMC from . .

Source codes are here. The details are here (in Japanese). statmodeling.hatenablog.com I'd like to supplement the answer to the question. Simulated annealing can reach optimum, but replica exchange MCMC can get the distribution. That is a big difference.

Michael Betancourt "Scalable Bayesian Inference with Hamiltonian Monte Carlo"

He started his talk from an overview of tall data and wide data, then stepped forward to Hamiltonian Monte Carlo (HMC) which can do efficient sampling of high dimensional models (i.e. many parameters). His introduction of HMC using gravity was really nice. Xiangze who was one of the participants posted a good blog about that (in Japanese). xiangze.hatenablog.com NUTS is an algorithm that can tune the parameters in HMC automatically. The details of NUTS seemed to be cut because it was hard for us to comprehend. He also introduced Adiabatic Monte Carlo (arXiv) maybe for my presentation. Thank you! The animation was so cool.

Michael Betancourt "Some Bayesian Modeling Techniques in Stan"

He talked about linear models, GLM and hierarchical models gradually and kindly. Centered vs. non-centered parameterization was very interesting. Centered models are ordinary hierarchical models, and non-centered models are so-called reparameterized (scale-separated) hierarchical models. Please refer to the above video or the JSS paper of Stan if you want to know details.

I'd like to quote part of JSS paper of Stan.

Betancourt and Girolami (2013) provide a detailed analysis of the benefits of centered vs. non-centered parameterizations of hierarchical models in Euclidean and Riemannian HMC, with the conclusion that non-centered are best when data are sparse and centered when the data strongly identifies coefficients.

However, it is generally difficult to measure the sparseness of data, so the recommendation could be the following (from the JSS paper).

We recommend starting with the more natural centered parameterization and moving to the non-centered parameterization if sampling of the coefficients does not mix well.

It seems that we can also use non-centered parameterization partially (i.e. a stratum or a group).

Thank you very much again, DWANGO, our stuffs, Hiroki and Michael Betancourt!