StatModeling Memorandum

StatModeling Memorandum

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

Stanの関数を使ってRを拡張して高速化する

C++に自動で変換される)Stanの関数を使ってRを拡張できる機能が、Stan/RStanの2.16で実装開始されて2.17でほぼ完成しました。Rを高速化するためにC++(とRcpp)はあまり書きたくないけれど、Stanの関数なら書いてもいいよという僕得な機能です。この記事ではその方法を簡単に紹介します。

元にした資料はRStanの開発者であるBenさんがStanCon2018で発表したこちらの資料です。

ここでは例として、以下の2つの関数をRで使えるようにしましょう。

  • 1) 機械学習分野でおなじみのlog_sum_exp関数
    • 引数はN個の正の実数
  • 2) データにemax modelという曲線をあてはめた場合の対数尤度を返す関数
    • 引数はデータ(N個のXYのペア)とパラメータの値

手順は簡単で以下だけです。

  • functionsブロックだけ書いたstanファイルを用意する
  • R側でrstan::expose_stan_functions()を呼び出す

具体的には、以下の内容を例えばmy_functions.stanという名前で保存します。

functions {
  real stan_log_sum_exp(vector x) {
    return(log_sum_exp(x));
  }

  real emax_model_ll(vector Y, vector X, real e0, real emax, real ed50, real sigma) {
    int N = num_elements(Y);
    vector[N] mu = e0 + emax*X ./ (ed50 + X);
    real ll = normal_lpdf(Y | mu, sigma);
    return(ll);
  }
}
  • 3行目:log_sum_exp関数の方はすでにStanに用意されているものを使います。
  • 6行目:接尾の_llはlog likelihoodの意味です。
  • 7行目:たしかStan 2.14ぐらいから変数の宣言と定義を同時にできるようになりました。ここではデータポイントの数NYから動的に取得しています。

次にR側でrstan::expose_stan_functions()を呼び出します。

rstan::expose_stan_functions('my_functions.stan')

するとfunctionsブロックに入っていた関数は全てR側で使えるようになります。

stan_log_sum_exp(c(0.1, 0.05, 0.8))
emax_model_ll(Y=c(1,10,12,24), X=c(0,5,10,10), e0=0.8, emax=30, ed50=5.1, sigma=3)

StanのコードはC++のコードに変換されるのでfor文はRと比べて非常に高速です。また、Stanのマニュアル(特にBuilt-In Functionsの章以降)を見ていただくとよいのですが、Stanはマニアックな分布の対数尤度を求める関数や、行列演算の関数(Eigen由来)や、log1m_expなどの指数対数を安定に計算する関数が豊富です。

これらのStanの長所を非常に簡単にRのコーディングに取り込めるのは最高です。まあ、少しばかりStanの関数やデータ型に慣れている必要がありますが。

Enjoy!

統計・数学・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]にならって区別しないでそう呼びます。