StatModeling Memorandum

StatModeling Memorandum

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

統計・数学・R関連で用途別のオススメ書籍

比較的読みやすい本を中心に紹介します。2016年も同様の記事を書きましたが新しい本も出たので更新しました。今後はこのページを更新します。

数学入門

かなり前に大学卒業した人が数学に再入門しようと思っても、いきなり大学数学の問題集は解けません。もしくは式変形を眺めて分かった気になるだけで成長がありません。そこで高校数学の「黄色チャート」で特訓です。独力で式変形する練習をします。例題と重要例題だけ解きまくるのがオススメです。I・II・IIIは続き物なので、Iが簡単すぎるならIIやIIIへ進むこともできます。Aは場合の数と数列、Bはベクトル、Cは行列と確率分布が重要です。勉強と思わず、電車内で遊ぶパズルと思って純粋に楽しんだらいいと思います。

改訂版チャート式解法と演習数学1+A

改訂版チャート式解法と演習数学1+A

チャート式解法と演習数学2+B

チャート式解法と演習数学2+B

チャート式 解法と演習数学3+C 改訂版

チャート式 解法と演習数学3+C 改訂版

微分積分

高校数学は完璧とすれば、大学の微積で必要なのは偏微分テイラー展開がしっかりできることと思ってます。そこで「新しい微積分<上>」をすすめます。

新しい微積分<上> (KS理工学専門書)

新しい微積分<上> (KS理工学専門書)

統計や機械学習の分野でよく使うけど、上記のような教科書にあまり登場しない話題として、ラグランジュの未定乗数法などがあります。そういう話題はネットや別の本で学びましょう。

線形代数

tensorflowなどのおかげで順伝播部分(行列積および行列とベクトルの積)さえ書ければ線形代数の知識はそこまでいらないんじゃないかという流れを感じます。しかし、主成分分析やトピックモデルなどの行列分解や、ガウス過程などのカーネル法のような様々なデータ解析の手法に一歩踏み込むと、きちんとした勉強が必要になります。理解しやすくて使いやすくて、統計や機械学習への応用を主眼においた線形代数の本はまだ見たことないです。機械学習シリーズとかで基礎から「The Matrix Cookbook」*1までを視野に入れた本が出てきてくれると良いのですが。

以下では個人的な好みの教科書を挙げておきます。「線形代数―基礎と応用」です。2x2行列の場合から徐々に説明していく点と、図形と関連づけて説明することが多い点が分かりやすくて好きです。

講座 数学の考え方〈3〉線形代数―基礎と応用

講座 数学の考え方〈3〉線形代数―基礎と応用

統計学入門

色々読んでみましたが、現在決定版と言えるものは存在しないように思えました。個人的には、シグマと積分の復習、場合の数・数え上げの方法、確率、確率変数、確率密度、度数分布とヒストグラム、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換、検定、散布図と箱ひげ図、回帰、相関あたりをRやPythonなどを使いながらシンプルに説明していく本があるといいと思うのですが、なかなかバランスのとれたいい本がありません。初歩の初歩しか説明してない、グラフが少ない、検定にページを割きすぎ、分厚い、ちょっと難しいなどの不満点があります。立ち読みして自分にあった本を選ぶのがいいと思います。

読み物では大村平先生のシリーズをおすすめします。最近、hoxo_m さんと僕がともに大村氏の本で確率統計をはじめたと分かって盛り上がりました。僕は「確率のはなし」「統計のはなし」「統計解析のはなし」「多変量解析のはなし」「実験計画と分散分析のはなし」「ORのはなし」を持っていますが、どれも読みやすくて面白いです。特に「統計解析のはなし」「実験計画と分散分析のはなし」「ORのはなし」は面白かった覚えがあります。

統計解析のはなし―データに語らせるテクニック (Best selected Business Books)

統計解析のはなし―データに語らせるテクニック (Best selected Business Books)

ORのはなし―意思決定のテクニック

ORのはなし―意思決定のテクニック

統計学演習

統計の実力向上のためには、「紙と鉛筆で手を動かして計算する」という経験がどこかで必要になります。そこで、シグマと積分はすでに慣れ親しんで、統計の初歩は少し知っている人に「統計学演習」をオススメします。

統計学演習

統計学演習

演習問題の難易度と内容が素晴らしいです。しかし、演習の前にある簡潔なまとめは説明がないので、復習として使うのがよいでしょう。

プログラミング言語としてのRに詳しくなりたい

「Rプログラミング本格入門: 達人データサイエンティストへの道」がオススメです。Rubyなど他のプログラミング言語に慣れているようなエンジニアに特に向いていると思います。前半はRの簡潔で良質なまとめ、後半はメタプログラミング・データベース操作・高速化手法など中級者向けの話題です。統計手法はあまり載っていません。僕はR言語徹底解説は8章から難しくなってきて読み切れなかったのですが、この本は読み切ることができました。

Rプログラミング本格入門: 達人データサイエンティストへの道

Rプログラミング本格入門: 達人データサイエンティストへの道

Rのデータ処理に詳しくなりたい

Rといえばデータ処理能力と作図能力がずば抜けています。データ処理については今のところ「Rではじめるデータサイエンス」をオススメします。翻訳者がRを普段から使っていないためなのか、日本語による記述はやや微妙な部分があると感じますが。

Rではじめるデータサイエンス

Rではじめるデータサイエンス

Rの可視化(ggplot2)を使ってみたい

今のところ「Rグラフィックスクックブック」がオススメです。カラーページが多くびびります。

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

プログラミング少しできる人が統計モデリングに詳しくなりたい

「StanとRでベイズ統計モデリング」をオススメします。StanEdwardPyroなどと同様の確率的プログラミング言語で、高次元のパラメータ空間からサンプリングを効率的に行えるのが特徴です。今後、確率的プログラミング言語で新しい独自モデルを試行錯誤していきたい人にとっても一読の価値があると思います。*2

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

Rを使って機械学習に詳しくなりたい

RやPythonを使って機械学習をやってみようという本は多いのですが、手法の背後にある考え方までもきちんと説明している本はあまり見たことありません。洋書になってしまいますが、「An Introduction to Statistical Learning with Applications in R」がオススメです。ここからpdfがダウンロードできます。著者の中はあのHastie先生とTibshirani先生もいます。この本をもとにしたオンライン授業もここから無料で受講することができます。阪大などの授業でも使われているようです。

因果推論や効果測定を詳しく知りたい

「岩波データサイエンス vol.3」がオススメです。統計を扱う上で非常に重要なのが因果関係ですが、この本を除いて読みやすい類書がほとんどありません。傾向スコアだけでなく、反実仮想を考えた回帰モデルやRCT(ランダム化比較試験)との対応がきちんと書いてあって勉強になりました。

岩波データサイエンス Vol.3

岩波データサイエンス Vol.3

Stanにもっと詳しくなりたい

Stanのマニュアルがオススメです。日本語化プロジェクトがあり、少し前のバージョンですが主要パートについては日本語訳もほとんど完成しています(ここ)。現在は原文に忠実にするよりも、より大胆にそぎ落として分かりやすい日本語にするのを目標にしてます。もともとの英語が難しく訳に苦戦している箇所もあるので、勉強がてら是非プロジェクトへの参加をお待ちしております。 github.com

*1:研究者にとってハンドブック的に有名な本。行列の微分とかが載ってます。pdf file

*2:Edwardはバグは少し多めで機能追加を優先する傾向があります。MCMCサンプリングは収束が悪かったり、バグがまだまだ見つかっており(これとか)で実用に耐えれる印象はありません。現在のところEdwardなら自動変分ベイズで推定するのがオススメです。

IGMRFの尤度におけるrankの減少分に関するメモ

以下の書籍を読んで、IGMRF(Intrinsic Gaussian Markov Ramdom Field)の尤度に関して自分の理解をまとめたメモです。

この発表資料の18ページにおいて、(観測モデル部分を除いた)IGMRFの対数尤度は以下に比例すると書きました(ただし \sigma_\mu \sigmaに、 \mu \boldsymbol{x}に読み替えてください)。

 

ここで nはnodeの数、 \textbf{Q} n \times nの精度行列*1でnodeのつながりの情報を反映していて、 \Deltaは線形制約に由来する \textbf{Q}のrankの減少分です。

1次元の1階階差のIGMRF

線状につながれたGMRFの場合、 \textbf{Q}は以下になります。

 

このように精度行列 \textbf{Q}が帯行列になってスパースになるところがGMRFの特徴です(分散共分散行列はスパースにならない)。

ここで、 \textbf{Q} \boldsymbol{1} = \boldsymbol{0}という線形の制約を満たすことに注意してください。 \boldsymbol{x}に定数を足して平行移動しても、 \boldsymbol{x}の要素の差しか尤度に関わってこないので、全体の尤度は不変であることがこの制約の由来です。数式で表現すると、 \boldsymbol{x} \boldsymbol{x} + \boldsymbol{1} \beta_0に平行移動しても、 (\boldsymbol{x} + \boldsymbol{1} \beta_0)^T \textbf{Q} (\boldsymbol{x} + \boldsymbol{1} \beta_0) \boldsymbol{x}^T \textbf{Q}\boldsymbol{x}と同じになります。 iを位置のインデックスとして x_i x_i + \beta_0に変えても同じという意味です。 \textbf{Q} \boldsymbol{1} = \boldsymbol{0}という制約のため、 \textbf{Q}のrankは1だけ減少して、 n-1となります。よって、IGMRFの対数尤度は以下になります。

 

1次元の2階階差のIGMRF

線状につながれたGMRFの場合、 \textbf{Q}は以下になります。

 

ここで、 \textbf{Q} \textbf{S}_{1} = \boldsymbol{0}という制約を満たします。  \textbf{S}_{1}は以下の行列です。

 

これは \boldsymbol{x}に直線を足しても差の差(例: (x_3 - x_2) - (x_2 - x_1))しか尤度に関わってこないので、直線の傾きがキャンセルされて全体の尤度は不変であることがこの制約の由来です。数式で表現すると、 \boldsymbol{x} \boldsymbol{x} + \textbf{S}_1 \boldsymbol{\beta}に変えても、 ( \boldsymbol{x} + \textbf{S}_1 \boldsymbol{\beta})^T \textbf{Q} ( \boldsymbol{x} + \textbf{S}_1 \boldsymbol{\beta} ) \boldsymbol{x}^T \textbf{Q} \boldsymbol{x}と同じになります。ここで \boldsymbol{\beta} = ( \beta_0 \ \ \beta_1 )^Tです。 iを位置のインデックスとして x_i x_i + \beta_0 + \beta_1 iに変えても同じという意味です。 \textbf{Q} \textbf{S}_{1} = \boldsymbol{0}という制約のため、 \textbf{Q}のrankは2だけ減少して n-2になります。よって、IGMRFの対数尤度は以下になります。

 

 d次元の k階階差のIGMRF

 d次元の k階階差のGMRFではrankはいくつになるでしょうか。結論を先に書くと、以下になります。

 

ここで、 nはnodeの数です。そこから二項係数を引いています。ここで、 \boldsymbol{x}をすべてのnodeを1列に並べたベクトルとします。すると、 \boldsymbol{x} \boldsymbol{x} + \textbf{S} \boldsymbol{\beta}に変えても尤度は不変となります。この \textbf{S}の列の数(= \boldsymbol{\beta}の要素数)がrankの減少分に相当します。これまでの例では \textbf{S}は、 \boldsymbol{1}だったり、 \textbf{S}_{1}だったりしました。これまでの例と同じように、 (x_{i_1}, x_{i_2}, ..., x_{i_n}) に加えても不変になる項を求めると以下になります。

 

ここから、 \boldsymbol{\beta}の要素数は次のように数えます。 \betaの添字の数が d個で、各添字の数値の範囲が 0 k-1で、添字の数値の合計は k-1以下となる制限があります。この状況は、部屋の仕切り d個と k-1個の玉の並び替えの場合の数と同じになるので、 \boldsymbol{\beta}の要素数は以下になります。

 

例えば2次元3階階差の場合、 (x_i, x_j) に加えても不変になる項は、以下になります。

  \beta_{00} + \beta_{10} i + \beta_{01} j + \frac{1}{2} \beta_{20} i^2 + \beta_{11} i j + \frac{1}{2} \beta_{02} j^2

 \boldsymbol{\beta}の要素数は6個になりますので、rankの減少分は6になります。例えば、 \beta_{10}は、以下の部屋の仕切りと玉の図に相当します。

 

Stanによる実装例

「StanとRでベイズ統計モデリング」の12.8節の例題(2次元1階階差の地図)のStanコードは以下になります。

data {
  int N;
  vector[N] Y;
  int I;
  int<lower=1, upper=N> From[I];
  int<lower=1, upper=N> To[I];
}

parameters {
  vector[N] r;
  real<lower=0> s_r;
  real<lower=0> s_Y;
}

model {
  target += - (N-1) * log(s_r) - 0.5 * dot_self(r[To] - r[From]) * inv_square(s_r);
  Y ~ normal(r, s_Y);
}

なお、先に \textbf{Q}を頑張って構築してquad_formなどを利用してモデルを書く方法は、上記のように書く方法より遅くなりました。ベクトル化がうまく効かないからだと思います。なお、12.8節のモデルと比べると \textbf{Q}の部分は等価ですが、外側のlog(s_r)に掛かる定数が異なります(12.8節のモデルがより \sigma_rが小さくなる方へペナルティが入っています)。

*1:厳密にはrankが落ちているので精度行列そのものとは言えませんが、ここでは[4]にならって区別しないでそう呼びます。

『ベイズ統計モデリング ―R,JAGS, Stanによるチュートリアル 原著第2版―』 John Kruschke著、前田和寛・小杉考司監訳

タイトルの本を頂きました。ありがとうございます。僕は原著を少し読んだことがあり、こちらで非常に評判が高い本です。翻訳にもかかわらず原著とほぼ同じ値段で購入できます。

先にJAGSになじみのない方へ説明しておきますと、JAGSはRコアメンバーの一人でもあるMartyn Plummer氏によってC++で開発されたMCMCソフトウェアです。Rから使うのが多数派ですが、PythonからもPyJAGSによって使うことができます。

複雑なモデルでなければStanより収束が早く、離散値をとるパラメータも使えるため、プログラミングがそんなに得意でない人がベイズ統計モデリングをはじめるには一番向いていると思います。最近、再び活発に開発され始めたようで、先日JAGS 4.3.0がリリースされました。

RからのJAGSの使い方としては、hikaru1122さんの記事がよくまとまっていてオススメです。 hikaru1122.hatenadiary.jp

特長

一言で言えば「ベイズ統計モデリングチュートリアル」としてオールインワンとなっている本です。

方向性が類似している和書は、豊田秀樹(2016)『はじめての統計データ分析』です。すなわち、伝統的な統計分析(t検定・分散分析、ロジスティック回帰などの一般化線形モデル)をすべてベイズ統計で焼きなおしていきましょうという方向性です。理論だけでなく、きちんとソースコード(R,JAGS,Stan)の解説もついています。さらに階層ベイズモデルやモデル比較についても章を割いて丁寧に説明しています。

しかし、そういう発展的な話題だけでもないのです。本の前半では、確率とは何か、条件付き確率・同時確率・独立、ベイズ統計の初歩についても説明しているほか、Rの基本、MCMCの基本(数式が多く割と難しめ章ですが竹林先生の訳がいい感じです)、BUGS/JAGSの基本(これまでの和書にはほとんどない情報!)、Stanの基本と、「チュートリアルとしてはこの本さえ読めば他は読まなくていい」というぐらいのオールインワンとなっています。そりゃそこまでやれば770ページになるし、鈍器にもなります。『はじめての統計データ分析』は比較的コンパクトにまとまっているために、初学者にとって本当に「はじめて」と言えるのか疑問な箇所がいくつかありましたが、この本はそのページ数を見ても分かる通り、手取り足取り非常に丁寧で「チュートリアル」のサブタイトルに偽りありません。

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

  • 作者: John K. Kruschke,前田和寛,小杉考司,井関龍太,井上和哉,鬼田崇作,紀ノ定保礼,国里愛彦,坂本次郎,杣取恵太,高田菜美,竹林由武,徳岡大,難波修史,西田若葉,平川真,福屋いずみ,武藤杏里,山根嵩史,横山仁
  • 出版社/メーカー: 共立出版
  • 発売日: 2017/07/22
  • メディア: 大型本
  • この商品を含むブログを見る

一点、類書と異なり注意を要する点は、この本ではベイズ信頼区間(信用区間, 確信区間)の定義として、両端からα/2%ずつ切って定義する方法ではなく、ほぼすべてで最高密度区間(Highest Density Interval; HDI*1)を使っている点です。こだわりと思います。HDIの方が都合の良い場合もあって、本には例として、峰が離れた二峰の事後分布になった場合に、区間が2つになって直感的に解釈できるとあって、なるほどと思いました。

また、この本は図もとても豊富で充実しています。特に以下のようなモデルの図が頻繁にあって分かりやすいです(この図は原著者の以下のブログ記事からの引用です http://doingbayesiandataanalysis.blogspot.jp/2015/12/lessons-from-bayesian-disease-diagnosis_27.html)。

書籍にはこのような可視化のコードも解説があります。しかし、ソースコードGitHubで管理されているわけではなく、以下のサイトからダウンロードするようです。

Software, with programs for book - Doing Bayesian Data Analysis

読者の方が非公式に更新しているリポジトリはあるようです。

GitHub - boboppie/kruschke-doing_bayesian_data_analysis: John K. Kruschke's Doing Bayesian Data Analysis: A Tutorial with R and BUGS

このあたり、訳書の方にも演習の答え(日本語)や図を含めたリポジトリがあるいいなぁと思いますが、どうでしょうか。期待しましょう。

数少ない不満点は、教養のある原著者で洋書にありがちな、流れを妨げる(理解に時間を要する)オヤジギャグ部分がたびたびあること。僕が訳者なら、原著者に「この部分、通じない可能性が高いのでカットしていいですか?」とか問い合わせてカットを試みたいです。まぁ、オヤジギャグが好きな読者もいるので、好みの問題かもしれません。

翻訳について

最近この手の技術関連の翻訳本を読む際に真っ先に調べることがあります。それは訳者の背景・経歴・アウトプットです。いわゆるプロ訳者の方でも、技術の鍛錬を怠っていると、日本語としては読みやすくとも不思議なことに中身がピンとこない翻訳になるのです。普段からその技術を追っていない(or使っていない)ために、「パッケージの方向性」、「著者の熱意がある部分」、「現実の問題への使われ方」などをよく理解できていないからだと思います。結果的にシャープでない曖昧な訳になります。その点、この本の訳者は数多くいますが、実際に自分たちの研究論文を書くために「ベイズ統計モデリング」を理解し使っている人々なので、曖昧な訳で逃げているような箇所はいまのところほとんど見当たらなかったです(オヤジギャグ部分はちょっとつらそうですが…)。また、訳語やトーンもほぼ一貫しており、監訳者二人の厳しい努力が垣間見えます。

*1:Highest Posterior Density Interval; HPDIという用語もあります