マルコフ連鎖モンテカルロ(MCMC)モデルのStanを使ったRとClojureによる解析

2022.10.17
機械学習技術 人工知能技術 デジタルトランスフォーメーション技術 確率的生成モデル 本ブログのナビ アルゴリズム 自然言語処理技術 深層学習技術 トピックモデル マルコフ連鎖モンテカルロ法 C/C++と機械学習 Clojureと関数プログラミング
マルコフ連鎖モンテカルロ(MCMC)モデルのStanを使ったRとClojureによる解析

ベイズ推定のような複雑な複率分布の積分計算を解析的に行うことは困難であり、マルコフ連鎖モンテカルロ(Markov Chain Monte Carlo;MCMC)法がよく用いられる。これは基本的な乱択アルゴリズムの一種で、最も単純なモンテカルロ法は、パラメータの候補として乱数を発生し、その乱数に対応する確率(積分)を計算するものだが、このような総当たり的な手法では、パラメータの数が増えると爆発的に計算量が増え効率的でないのに対して、完全にランダムではなく、ひとつ前の乱数を元に少しずつ確率の大きな値を探索(マルコフ確率場の探索)しつつ次の乱数次の乱数を発生するアルゴリズムとなる。詳しくは本ブログのマルコフ連鎖モンテカルロ(Markov Chain Monte Carlo;MCMC)法を参照のこと。

このMCMCを計算するためのフレームワークとして著名なものにSTANがある。STANはアンドリュー・ゲルマンらによってC++で開発され2012年に初版が提供されたものとなる。このStanに関してはR言語との統合であるRStan、Pythonとの統合であるPyStan、コマンドライン実行可能ファイルであるCmdStan、Clojureとの統合であるclj-stan等のさまざまな言語でのインターフェースを持つ。

以前、Rstanについては紹介したが、PythonでのMCMCの実装としてはPystan、PyMC3がある。また、Rstanは2020年から更新されておらず、StanチームもCmdStanを奨励しているらしいことから、今回はCmdStanのrのラッパーであるCmdStanrと同様にCmdStanのシンプルなラッパーであるclj-stanについて述べてみたいと思う。

Rを使ったCmdStanの実行

まずCmdStanrについて。Rのツールを起動させてCmdStanrをインストールする。

>install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

さらにその他の必要ライブラリをインストールする。

>install.packages("posterior")
>install.packages("bayesplot")

ダウンロードしたライブラリを読み込む

>library(cmdstanr)
>check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
>library(posterior)
>library(bayesplot)
>color_scheme_set("brightblue")

CmdStanをインストールする。

>install_cmdstan(cores=2)

ダウンロードが完了したら、パスとバージョンを確認してみる。

>cmdstan_path()
>cmdstan_version()

両方とも無事に出て来ればインストールは完了となる。

次に、Stanライブラリに入っているサンプルモデルを使ってみる。

>file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan")
>mod <- cmdstan_model(file)

内容は$print()メソッドで確認できる。

>mod$print()
data {
  int<lower=0> N;
  array[N] int<lower=0,upper=1> y; // or int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);  // uniform prior on interval 0,1
  y ~ bernoulli(theta);
}

上記のモデルでのdata、parametert、modelに関しては、RでのSTANを使ったベイズ推定を参照のこと。

早速MCMCを動作させてみる。

>data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))

>fit <- mod$sample(
  data = data_list, 
  seed = 123, 
  chains = 4, 
  parallel_chains = 4,
  refresh = 500 # print update every 500 iters
)

モデルもシンプルで、データ点数も少ないため、即座に計算は終了する。

Running MCMC with 4 parallel chains...

Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 1 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 1 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 2 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 2 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 3 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 3 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 4 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 4 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Chain 3 finished in 0.0 seconds.
Chain 4 finished in 0.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.3 seconds.

結果をサマライズするには、posteriorパッケージのsummarise_draws()を呼び出したsummary()メソッドを使う。これは最初の引数で要約する変数を指定し、それ以降の引数はposterior::summarise_draws()に渡され、どの要約を計算するか、マルチコアを使うか、などを指定することができる。

>fit$summary()
# A tibble: 2 × 10
  variable   mean median    sd   mad      q5    q95  rhat ess_bulk ess_tail
                         
1 lp__     -7.27  -7.00  0.709 0.344 -8.70   -6.75   1.00    1852.    2114.
2 theta     0.247  0.232 0.119 0.123  0.0804  0.466  1.00    1611.    1678.
Clojureでの実行

前述したClojisrを利用する。環境設定までは”ClojureとRの連携による統計的学習”を参照のこと。

project.cljファイルは以下の通り。

(defproject clj-stan "0.1.0-SNAPSHOT"
  :description "FIXME: write description"
  :url "http://example.com/FIXME"
  :license {:name "EPL-2.0 OR GPL-2.0-or-later WITH Classpath-exception-2.0"
            :url "https://www.eclipse.org/legal/epl-2.0/"}
  :dependencies [[org.clojure/clojure "1.11.1"]
                 [scicloj/clojisr "1.0.0-BETA20"]]
  :repl-options {:init-ns clj-stan.core}
  :jvm-opts ["-Dclojure.tools.logging.factory=clojure.tools.logging.impl/jul-factory"])

src/core.cljファイルは以下のようになる。

(ns clj-stan.core)

(require '[clojisr.v1.r :as r :refer
            [r eval-r->java r->java java->r java->clj java->native-clj
             clj->java r->clj clj->r ->code r+ colon require-r]]
          '[clojisr.v1.robject :as robject]
          '[clojisr.v1.session :as session])

(require-r '[cmdstanr]
           '[posterior]
           '[bayesplot])

(r "file <- file.path(cmdstan_path(), 'examples', 'bernoulli', 'bernoulli.stan')")
(r "mod <- cmdstan_model(file)")

(r "mod$print()")

(r "data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))")
(r "fit <- mod$sample( data = data_list, seed = 123, chains = 4, parallel_chains = 4, refresh = 500 )")

(println (r "fit$summary()"))

ライブラリの導入の部分がlabray(cmdstanr)→(require-r ‘[cmdstanr])となる以外は全てr関数でRのコードをそのまま評価することで、同じ動作を行えることが確認できる。

以前述べた”ClojureとPythonの連携と機械学習”で既存のpythonライブラリをClojure上で利用することもできるので、これらを使ってR⇄Clojure⇄Pythonの豊富なライブラリ処理をClojure上で組み合わせて全て利用することが可能となる。(例→Clojure(IOT等のリアルタイムデータ入力処理)→R(CmdStan)→Python(Pandasのライブラリで可視化)等)

コメント

  1. 確率関数の積分等に用いられるマルコフ連鎖モンテカルロ | Deus Ex Machina より:
    2023年6月6日 2:13 PM

    […] マルコフ連鎖モンテカルロ(MCMC)モデルのStanを使ったRとClojureによる解析 […]

    返信
  2. ベイズ深層学習の概要と適用事例及び実装例 | Deus Ex Machina より:
    2023年8月31日 4:46 AM

    […] 確率的プログラミング: “Clojureを用いた確率的プログラミング(Probabilistic Programming)“でも述べている確率的プログラミングは、確率モデルを記述するための高水準のプログラミング言語を使用してモデルを定義し、モデルの事後分布を推定する手法となる。”マルコフ連鎖モンテカルロ(MCMC)モデルのStanを使ったRとClojureによる解析“にも述べているStanやPyroなどのライブラリが確率的プログラミングのために利用される。 […]

    返信
タイトルとURLをコピーしました