StatModeling Memorandum

StatModeling Memorandum

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

スパースモデルではshrinkage factorの分布を考慮しよう ~馬蹄事前分布(horseshoe prior)の紹介~

ベイズ統計の枠組みにおいて、回帰係数の事前分布に二重指数分布(ラプラス分布)を設定し回帰を実行してMAP推定値を求めると、lassoに対応した結果になります。また、回帰係数にt分布を設定する手法もあります。これらの手法は「shrinkage factorの分布」という観点から見ると見通しがよいです。さらに、その観点から見ると、馬蹄事前分布が魅力的な性質を持っていることが分かります。この記事ではそれらを簡単に説明します。

なお、lassoそのものに関しては触れません。岩波DS5がlassoを中心にスパースモデリングを多角的に捉えた良い書籍になっているので、ぜひそちらを参照してください。

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

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

  • 発売日: 2017/02/16
  • メディア: 単行本(ソフトカバー)

参考文献

  • [1] C. Carvalho et al. (2008). The Horseshoe Estimator for Sparse Signals. Discussion Paper 2008-31. Duke University Department of Statistical Science. (pdf file)
  • [2] J. Piironen and A. Vehtari (2016). On the Hyperprior Choice for the Global Shrinkage Parameter in the Horseshoe Prior. arXiv:1610.05559. (url)
  • [3] C. Carvalho et al. (2009). Handling Sparsity via the Horseshoe. Journal of Machine Learning Research - Proceedings Track. 2009;5:73–80. (pdf file)
  • [4] Epistemology of the corral: regression and variable selection with Stan and the Horseshoe prior:PyStanの開発者による、PyStanでの実装例とlassoとの比較です。

horseshoe+ priorなんてのも提案されています。

理論編

分割された正規分布

まず分割された正規分布の復習から。例えば、PRMLの2.3.2項に記述があります。以下の多変量正規分布があるとします。

f:id:StatModeling:20201106171132p:plain

ここで、以下のように2つに分割します。

f:id:StatModeling:20201106171136p:plain

f:id:StatModeling:20201106171141p:plain

f:id:StatModeling:20201106171145p:plain

すると条件付き分布は以下になります。

f:id:StatModeling:20201106171148p:plain

f:id:StatModeling:20201106171153p:plain

shrinkage factorの導出

上記の[2]を参考にしました。丁寧な導出が見つからなかったので、以下は手計算で求めました。間違っていましたら連絡ください。

以下の線形モデルを考えます。

f:id:StatModeling:20201106171157p:plain

ここで \overrightarrow{\beta}が回帰係数のベクトルで \bf{X}が説明変数の行列です。多変量正規分布で表現すると以下になります。

f:id:StatModeling:20201106171202p:plain

ここで、 \bf{\Lambda}は以下の対角行列です。

f:id:StatModeling:20201106171205p:plain

 * * *

さて、ここで以下の \overrightarrow{z}を考えます。

f:id:StatModeling:20201106171209p:plain

同時分布の対数は以下になります。

f:id:StatModeling:20201106171213p:plain

これは \overrightarrow{z}の要素の2次関数なので、 p(\overrightarrow{z})正規分布になります。PRMLの2.3.3項とほぼ同じように、2次の項と1次の項にわけて考えることで、その正規分布の精度行列と平均ベクトルを求めることができます。まずは2次の項を整理すると、

f:id:StatModeling:20201106171217p:plain

f:id:StatModeling:20201106171221p:plain

となり、精度行列 \bf{T}が求まります。1次の項はないので、平均ベクトルは以下になります。

f:id:StatModeling:20201106171226p:plain

以上より、前述の「分割された正規分布」の公式より、 \overrightarrow{y}が与えられたもとでの \overrightarrow{\beta}の分布は以下になります。

f:id:StatModeling:20201106171229p:plain

ここで、平均ベクトルは以下のように計算できます。

f:id:StatModeling:20201106171234p:plain

f:id:StatModeling:20201106171239p:plain

f:id:StatModeling:20201106171242p:plain

f:id:StatModeling:20201106171246p:plain

ここで、 \overrightarrow{\hat{\beta}} \overrightarrow{\beta}の分布を考えない場合の最尤推定の解で以下です。

f:id:StatModeling:20201106171249p:plain

 * * *

さて、ここで各説明変数列が相関がなく、平均がゼロで、分散が1とします。すなわち、以下を仮定します。

f:id:StatModeling:20201106171252p:plain

すると、平均ベクトルは以下のように変形できます。

f:id:StatModeling:20201106171256p:plain

f:id:StatModeling:20201106171300p:plain

対角行列になるので要素ごと表すと、以下になります。

f:id:StatModeling:20201106171304p:plain

ここで、 \kappa_jは以下となり、shrinkage factorと呼ばれます。

f:id:StatModeling:20201106171307p:plain

shrinkage factorが0に近いと係数が最尤推定値に近くなり(shrinkageしてない)、1に近いと0に近くなります(shrinkageする)。

馬蹄事前分布(horseshoe prior)

上のモデルで各 \lambda_jがある確率分布に従うと仮定します。この確率分布によって、元の \beta_jにどんな確率分布を設定したのと等価になるのか、そしてshrinkage factorの分布がどのようになるのか手計算で求めることができます(参考: 確率変数の変数変換とヤコビアンに慣れる - StatModeling Memorandum)。以下に結果だけ対応表として載せておきます(手計算で求めましたが間違っていたらすみません)。

法名  \beta_jの分布  \lambda_jの分布  \kappa_jの分布
lasso double-exponential f:id:StatModeling:20201106171311p:plain f:id:StatModeling:20201106171322p:plain
- Student-t*1 f:id:StatModeling:20201106171315p:plain f:id:StatModeling:20201106171326p:plain
horseshoe f:id:StatModeling:20201106171334p:plain f:id:StatModeling:20201106171319p:plain *2 f:id:StatModeling:20201106171330p:plain

shrinkage factor( \kappa_j)の分布を描くと以下になります。

f:id:StatModeling:20201106171128p:plain

horseshoeの場合の \kappa_jの分布形が馬蹄に似ているので馬蹄事前分布(horseshoe prior)と名づけられました。なお、 n \sigma^{-2} \tau^{2} = 1のときはベータ分布 Beta(0.5, 0.5) に一致します。

馬蹄事前分布のようにshrinkage factorが0か1を生成しやすいということは、「0につぶしやすい一方で、0につぶれない係数は0に近づくような振る舞いをしないで最尤推定値に近づく」ことになります。係数にメリハリがあるわけです。これに対し、lassoの場合は n \sigma^{-2} \tau^{2}が1もしくは10だとあまり0につぶす傾向はなく、またつぶれなかった係数は少しだけ最尤推定値に近づく感じになります。 n \sigma^{-2} \tau^{2} = 0.1だと0につぶす傾向はあるのですが、つぶれなかった係数も0に近づけてしまいます。

この結果、C. Carvalhoらはシミュレーション実験から「一部だけ影響のない説明変数があるようなシミュレーションデータに対しては馬蹄事前分布の方が予測の性能がよい」と主張しています([3])。さらに馬蹄事前分布は二重指数分布(ラプラス分布)とは異なり至るところで微分可能なので、偏微分を使うStanでの推定が安定であることから、Stanの開発メンバーはスパースモデリングには二重指数分布よりも馬蹄事前分布の使用をすすめています。

Stanによる実装例

[4]の実装に尽きています。すなわち以下です。簡単です。

data {
  int<lower=0> N;
  int<lower=0> D;
  matrix[N,D] X;
  vector[N] Y;
}

parameters {
  vector[D] beta;
  vector<lower=0>[D] lambda;
  real<lower=0> tau;
  real<lower=0> sigma;
}

model {
  lambda ~ cauchy(0, 1);
  tau ~ cauchy(0, 1);
  for (d in 1:D)
    beta[d] ~ normal(0, lambda[d] * tau);
  Y ~ normal(X * beta, sigma);
}

Enjoy!

*1:参考: http://stats.stackexchange.com/questions/52906/student-t-as-mixture-of-gaussian

*2:正の部分だけ定義された半コーシー分布

『Pythonで体験するベイズ推論 ―PyMCによるMCMC入門―』の書評

特長

Pythonユーザが待ちに待ったPythonによるMCMC本ではないでしょうか。原著タイトルが『Bayesian Methods for Hackers』だけあって、プログラマ・エンジニア向きだと思います。数式はびっくりするほど出てこない代わりに、Pythonコードは非常にたくさんでてきます。そしてPyMCの使い方が基礎から説明してあって丁寧です。自分でコーディングする際は原著のGitHubリポジトリを活用しましょう(なんとStarが10000個を超えてる!)。

Pythonで体験するベイズ推論 PyMCによるMCMC入門

Pythonで体験するベイズ推論 PyMCによるMCMC入門

  • 作者: キャメロン・デビッドソン=ピロン,玉木徹
  • 出版社/メーカー: 森北出版
  • 発売日: 2017/04/06
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログを見る
購入を迷っている人の一番の心配は、本書のPyMCのバージョンが1つ前のPyMC2であることだと思います。しかし、そこは気にせず、まずは安定でドキュメントも多いPyMC2で勉強したらいいと思います。なぜなら、MCMCを使ったモデリングに慣れていない人のボトルネックは、現象を数式やコードに落とすモデル化という営みそのものに慣れていないことだからです。モデル化に慣れてPyMC2に慣れたら、PyMC3に移行するのはエンジニアの方なら大して時間かからないでしょう。僕自身、BUGS言語(WinBUGS, JAGS)に慣れてからStanに移行しましたが、おおよそ1週間で移行できました。迷わずPyMC2でモデリングライフを開始したらいいと思います。なお、原著のGitHubリポジトリにはPyMC3のコードも含まれています。

ちなみに2015年10月刊行の岩波データサイエンス Vol.1では、渡辺さんの頑張りのおかげでPyMC3の解説(20ページ程度)になっております。ご参考までに。

本書で最高だった点

5章です。秀逸すぎる。この章では事後分布と損失関数を組み合わせて意思決定する、いわゆるベイズ決定を行います。著者の金融分野における経験が生かされていて、さまざまな損失関数が出てきて有益情報がいっぱいです。以下では例を挙げますが、この他にも多数あります。詳しくは本を読んでください。

政治評論家向きの損失関数

 L(\theta, \hat{\theta})=\frac{|\theta - \hat{\theta}|}{\theta (1 - \theta)}

ここで、 \theta \hat{\theta}はともに [0,1]の範囲です。真の値 \thetaが0,1に近い場合には、優柔不断(ハッキリ0/1で答えない)だと損失が大きくなるような損失関数になっているとのことです。

株価の予測における下振れリスクを表す損失関数

下振れリスクとは「符号の異なる間違った方向に大きく予測すること」で、上振れリスクは「正しい方向に大きく予測すること」。下振れリスクの方が避けたいので、符号をまたいだところから損失がぐっと増えるような非対称な関数形になっているとのことです。図5.5(p.158)参照。

その他にも良かった点

  • ベイズの定理や大数の法則などの基礎となる理論も独特の切り口で扱っていて勉強になりました。
  • 頻繁に出てくる例題が軽妙でよく練られていて面白いです。例を挙げます。
    • ベイズの定理の説明で、農村における性格がこまかい人が将来司書になるか農家になるか問題。
    • 導入の1章において、著者に送られてきたメッセージの数の日次データに対する変化点検出。
    • Kaggleのダークマターハローの座標を予測する問題。銀河のプロットも面白い!
  • 事前分布の選び方で、統計の知識が全くない人から信念を聞き出して分布にする方法(ルーレット法)。

少し残念だった点

  • これはPythonのデータ解析・機械学習の本にありがちですが、可視化のコードの分量がやや多い気がします。matplotlibのせいかもしれません。仕方ないのかもしれません。
  • ベイズ推定の本領発揮とも言える、階層モデルなどのやや高度なモデルは扱っていないのは少し残念です。
  • MCMCの収束に関してやや議論が甘いと思います。他のMCMCソフトウェアの収束判定では \hat{R}(アールハット)がよく使われており、かなりよい指標だと思いますが、PyMCはなぜかchainを複数流して \hat{R}を求めることはしません。そして、書籍からもそこは抜けており、いくつかの解析で収束が怪しそうなものがあります。

{Rcpp}と{RcppGSL}を活用した数値積分の例

GSLはGNU Scientific Libraryの略で, 広い用途の数値計算向けのC言語のライブラリです。20年前から開発されており、まだリリースされ続けています。個人的な印象としては、めちゃくちゃチューニングされてはいないが、長年の開発のおかげで安定しており、マニュアル(htmlはこちら、pdfはこちら)も充実していて使いやすいイメージです。

この記事では{Rcpp}{RcppGSL}パッケージを通してRからGSLを使う例(数値積分)を挙げます。このような{Rcpp}を使って数値計算を行う類似パッケージには{RcppNumerical}が挙げられます。

なお、RからC++を使うためのパッケージである{Rcpp}自体の解説はしません。{Rcpp}の使い方とメリットに関しては、最近非常によい本が出版されました。以下の本の24章を参照するとよいと思います。

インストール

GSLのインストール

Linuxの場合

僕の動作環境はRは3.3.1、gccは6.1.1、gslは2.3です。

パッケージマネージャのようなものからインストールすることもできますが、僕は諸事情でソースからインストールしました。Linuxで典型的な手順を踏めばインストールできます。例えば、公式ページからtar.gzをダウンロードして、解凍してconfiguremakemake installしてしばらく待ちます。

その後、.bashrcなどで「gslをインストールしたディレクトリ」を以下のように環境変数LD_LIBRARY_PATHに追加しておきます。

export LD_LIBRARY_PATH=".:「gslをインストールしたディレクトリ」/lib:$LD_LIBRARY_PATH"

Macの場合

持ってないので分かりません。すみません。

Windowsの場合

僕の動作環境はRは3.3.1、RtoolsはRtools33またはRtools34です。

最新版のGSLでなくてもよければ、Ripley先生が提供しているサイトから、Windows用にprebuildされたGSL(local323.zip)をダウンロードして、、、と思ったら、先週ぐらいから403 Forbiddenになっているようです。終了です。*1Ripley先生が退職されたあとも今まで残っていたのが奇跡なのかもしれません。Cygwinからだとインストールできるかもですが、Rtoolsとの連携もしんどそうで頑張る元気がありません。素直にLinuxMacを使った方がよさそうです。

{Rcpp}と{RcppGSL}のインストール

いつも通りでOKです。すなわち、Rを起動してinstall.packages(c('Rcpp', 'RcppGSL'))でいけます。

テスト

テスト用のコードは作者のページにあるものを使います。C++ファイルとしてtest0.cppというファイル名で以下の内容で保存します。[[Rcpp::depends(RcppGSL)]]に注意です。

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

// [[Rcpp::export]]
Rcpp::NumericVector colNorm(const RcppGSL::Matrix & G) {
  int k = G.ncol();
  Rcpp::NumericVector n(k);
  for (int j = 0; j < k; j++) {
    RcppGSL::VectorView colview = gsl_matrix_const_column(G, j);
    n[j] = gsl_blas_dnrm2(colview);
  }
  return n;
}

さらに、別途Rファイルとしてtest0.Rというファイル名で以下の内容で保存します。

library(Rcpp)
library(RcppGSL)

sourceCpp('test0.cpp')

set.seed(123)
m <- matrix(rnorm(4), ncol=2)
print(colNorm(m)) #=> 0.6059 1.5603
  • 4行目:ここで先ほど作成したC++ファイルへのパスを記入します。同じディレクトリに入っていれば上のままでOKです。

Rコンソールを起動してそのディレクトリに移動し、source('test0.R')を実行して無事に0.6059 1.5603が出力されれば完了です(小数点何位まで出力されるかはRのオプションの指定によって異なると思います)。

RcppGSL活用例:数値積分

Rに数値積分するintegrate関数があるのに、わざわざ自分で数値積分を書くバカがどこにおるの?と思われるかもしれませんが、近頃おじさんはRのintegrate関数が機能しない場面にしょっちゅう出くわします。例を挙げます。

 \int_{-\infty}^{\infty} (0.2 + min(x, -0.3-x)) \times Normal(x-\mu, \sigma) dx

ここで、 \mu \sigmaは値が与えられています。例えば、 \sigma=0.24とし、 \mu -3から 3の間のいくつか選んで算出してみます。数値積分の範囲は {-\infty}から {\infty}まできちんと計算しなくても、 \mu - 6 \sigmaから \mu + 6 \sigmaまでの積分でほぼ十分とします。

Rだけで算出する場合は以下になります。

sd <- 0.24
integrand <- function(x, mu) (0.2 + min(x, -0.3-x)) * dnorm(x, mean=mu, sd=sd)
g1 <- function(mu) integrate(integrand, mu-6*sd, mu+6*sd, mu)$value
g2 <- function(mu) pracma::romberg(integrand, a=mu-6*sd, b=mu+6*sd, mu=mu)$value

x <- seq(-3, 3, len=101)
m1 <- sapply(x, g1)
m2 <- sapply(x, g2)

g1integrate関数で、g2は比較ためのRomberg法で数値積分した結果です。{pracma}パッケージのromberg関数で実行できます。この結果を可視化すると以下になります。

f:id:StatModeling:20201106175014p:plain

値がだいぶずれていますが、romberg関数の方が正しい値です。それならばromberg関数をいつも使えばいいじゃんと思うかもしれませんが、計算に時間がかかるのと、ここでは挙げませんがromberg関数だとかえって数値計算がうまくいかない場合もあります。そこで、{RcppGSL}パッケージからで数値積分します。

テストの場合と同様に、C++ファイルとしてgsl_cquad.cppというファイル名で以下の内容で保存します。

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_randist.h>

struct my_f_params { double mu; double sd; };

double integrand(double x, void *p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double mu = params->mu;
  double sd = params->sd;
  double res = (0.2 + fmin(x, -0.3-x)) * gsl_ran_gaussian_pdf(x-mu, sd);
  return(res);
}

// [[Rcpp::export]]
Rcpp::NumericVector integrate_my_cquad(Rcpp::NumericVector xVec, double sd) {
  Rcpp::NumericVector resVec(xVec.size());
  gsl_integration_cquad_workspace *work = gsl_integration_cquad_workspace_alloc(1e2);

  double result(0.0);
  double abserr(0.0);
  size_t nevals(1e2);
  double mu(xVec[0]);

  gsl_function F;
  F.function = &integrand;
  struct my_f_params params = { mu, sd };

  for (int i=0; i < xVec.size(); ++i) {
    mu = xVec[i];
    params.mu = mu;
    F.params = &params;
    int res = gsl_integration_cquad(&F, mu-6*sd, mu+6*sd, 1e-8, 1e-6, work, &result, &abserr, &nevals);
    resVec[i] = result;
  }

  gsl_integration_cquad_workspace_free(work);
  return resVec;
}
  • 3行目: GSLの数値積分の機能を使用する場合に必要となるヘッダーファイルです。
  • 4行目: GSLの確率分布を使用する場合に必要となるヘッダーファイルです。
  • 8~14行目: 被積分関数です。x積分変数で、pはその他の変数(パラメータ)です。RcppGSLで数値積分する場合、被積分関数はR側で作って渡すのではなく、C++側で定義しないと高速にならないので注意です。なお、gsl_ran_gaussian_pdf正規分布確率密度関数です。
  • 6行目: いま説明したpに値を渡すための構造体です。
  • 17~40行目: 数値積分を実行する自作関数を定義しています。RのvectorxVecに渡され、そのxVecの各値ごとに数値積分を繰り返し実行して、その結果をRのvectorとなるresVecで返します。数値積分にもいろいろ手法があります。特異点がある場合にRのintegrate関数がこける印象があることから、QAGSCQUADを検討した結果、CQUADの方が安定で速いケースが多かったので、この記事ではCQUADを使います。
  • 34行目: gsl_integration_cquadは数値積分をCQUADで実行する関数です。引数の説明についてはマニュアルを参照してください。引数の一つであるworkは作業領域になります。19行目のgsl_integration_cquad_workspace関数でその領域を確保し、38行目のgsl_integration_cquad_workspace_free関数でその領域を開放します。

C++ファイルで定義したintegrate_my_cquadを使うRスクリプトは以下になります。

library(Rcpp)
library(RcppGSL)

sourceCpp('gsl_cquad.cpp')

sd <- 0.24
x <- seq(-3, 3, len=101)
m3 <- integrate_my_cquad(x, sd)

計算の実行速度は非常に速く、また安定で、結果も上記のRomberg法とほぼ一致します。僕の実行環境で1000回ほどintegrate_my_cquad関数と同じ処理をして比較したところ、GSL版はRのromberg関数の約150倍、Rのintegrate関数の約1.5倍高速でした。

数値積分の簡単な使い方はこのサイトを参考にしました。

*1:portingのページはかろうじて残っていますが、、