StatModeling Memorandum

StatModeling Memorandum

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

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

R Stan

（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個のXとYのペア）とパラメータの値

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

  • 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ぐらいから変数の宣言と定義を同時にできるようになりました。ここではデータポイントの数NをYから動的に取得しています。

次に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!

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

Stan 理論

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

[asin:B00YBV6YLI:detail]

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

f:id:StatModeling:20201106165153p:plain

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

1次元の1階階差のIGMRF

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

f:id:StatModeling:20201106165157p:plain

このように精度行列 \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の対数尤度は以下になります。

f:id:StatModeling:20201106165201p:plain

1次元の2階階差のIGMRF

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

f:id:StatModeling:20201106165204p:plain

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

f:id:StatModeling:20201106165209p:plain

これは \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の対数尤度は以下になります。

f:id:StatModeling:20201106165212p:plain

d次元の k階階差のIGMRF

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

f:id:StatModeling:20201106165216p:plain

ここで、 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}) に加えても不変になる項を求めると以下になります。

f:id:StatModeling:20201106165223p:plainf:id:StatModeling:20201106165220p:plain

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

f:id:StatModeling:20201106165223p:plain

例えば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}は、以下の部屋の仕切りと玉の図に相当します。

f:id:StatModeling:20201106165227p:plain

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著、前田和寛・小杉考司監訳

書評 BUGS

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

先に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版

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

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

f:id:StatModeling:20201106165942p:plain

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

Doing Bayesian Data Analysis - Software, with programs for book

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

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という用語もあります