gretaパッケージを活用したベイズモデリングの紹介
概要
今回はRのgretaパッケージを使ったベイズモデリングを紹介します。
公式サイトは以下。
https://greta-dev.github.io/greta/index.html
公式のサンプルコードは以下にのっています。
gretaとは
gretaパッケージはどんなパッケージかと言うと、tendorflowをバックで動かしてモデリング可能なパッケージです。
色々出来るようなのですが、今回は公式サイトでも紹介されているようなMCMCを用いたベイズモデリングを試します。
特徴としては以下の3点があげられていて、
simple
- Rで書かれていて、BUGSやStanを覚える必要がない
scalable
extensive
- R関数を簡単にかける、また他のRパッケージを活用できる
ということらしいです。
これだけ読むと確率的プログラミング言語として、BUGSやStanが意識されている気がします。
実際、サイトのexampleを読むと、BUGSとStanのコードと比較した物が乗っており、gretaでは短くシンプルに書くことが出来ることを示されています。
インストール方法
まずは以下記事を参考にgretaパッケージを入れます。
install.packages("greta")
次にtensorFlowが必要なので、インストールしましょう。
インストール用の関数はgretaパッケージを読み込むと用意されていますので、以下を実行。
library(greta) install_tensorflow()
また以下のようなお洒落な図を書くためにDiagrammeRを入れます(図の意味は後述)。
install.packages('DiagrammeR')
これでgretaを使う準備ができました。
モデリング
今回はモデリングではirisデータを使って二つのモデルを試したいと思います。
具体的には、
Sepal.Lengthの平均とばらつきを種ごとに推定するモデル
Sepal.Lengthを従属変数、Petal.Lengthを説明変数、そして切片と傾きが種ごとに異なるランダム切片・係数モデル
です。
1. Sepal.Lengthの平均とばらつきを種ごとに推定するモデル
まずは種ごとのSepal.Lengthの平均とばらつきを推定するシンプルなモデルを試してみます。
以下コードを実行すると、
# library library(greta) library(tidyverse) # data x <- as_data(iris$Petal.Length) y <- as_data(iris$Sepal.Length) Species <-as.numeric(iris$Species) n_species <- length(unique(Species)) # variables and priors mu <- normal(0, 10, dim = n_species) sd <- uniform(0, 10, dim = n_species) coef <- normal(0, 10, dim = n_species) # operations mean <- mu[Species] + coef[Species] * x sd_species <- sd[Species] # likelihood distribution(y) <- normal(mean, sd_species) # defining the model m <- model(mu, sd, coef) # plotting plot(m)
こういった図が表示されると思います。
この図の各図形の意味は公式サイトにも乗っていますが、
白抜きの四角：データ
灰色の丸：変数
ひし形：確率分布
小さい丸：演算
線：決定した値
点線：確率値
となります。
先ほどプロットした図では種ごとに異なるといった階層性がうまく表現されていないので少し分かりづらいですが、 実際には種ごとに平均とばらつきが異なるといったモデリングしています。
図は階層性が表現されていないことを目をつぶるととてもわかり易いですね。
数式で表現すると、
なります。
ここでkは種の数なのでk=3です。
ここまででモデルが表現できたので、サンプリングをしてみましょう。
コードは簡単で、
# sampling draws <- mcmc(m, n_samples = 1000) summary(draws)
を実行するだけです。
結果は、
Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu1 5.0000 0.05439 0.001720 0.003638
mu2 5.9407 0.07065 0.002234 0.003926
mu3 6.5972 0.09385 0.002968 0.005132
sd1 0.3602 0.03853 0.001218 0.003890
sd2 0.5246 0.05261 0.001664 0.003478
sd3 0.6497 0.06765 0.002139 0.003892
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
mu1 4.8860 4.9583 5.0055 5.0391 5.0930
mu2 5.7902 5.8944 5.9402 5.9927 6.0839
mu3 6.4143 6.5404 6.5911 6.6624 6.7751
sd1 0.2980 0.3311 0.3563 0.3860 0.4508
sd2 0.4460 0.4838 0.5189 0.5601 0.6306
sd3 0.5319 0.5980 0.6436 0.6903 0.7884
とてもわかり易い結果ですね。
今の所Rhat値のような収束したかどうかの値を見る方法が分かっていないので（見れないのでしょうか）、 収束判断が難しいですが、 以下コードでサンプリングの結果や、信用区間を可視化することが出来ます。
library (bayesplot) #result plot mcmc_trace(draws) mcmc_intervals(draws)
ちなみに、mcmc関数は引数が以下のようになっているので、サンプリング手法等は変えることが出来るようです。
mcmc(model, method = c("hmc"), n_samples = 1000, thin = 1, warmup = 100,
verbose = TRUE, pb_update = 10, control = list(),
initial_values = NULL)
2. Sepal.Lengthを従属変数、Petal.Lengthを説明変数、そして切片と傾きが種ごとに異なるランダム切片・係数モデル
ランダム切片・係数モデルに関しては、コードのみを貼っておきますが、ほとんど先ほどと変わらず、 Petal.Lengthを説明変数として使って、切片と傾きが種ごとで異なるようにモデリングをしています。
# library library(greta) library (bayesplot) library(tidyverse) # data x <- as_data(iris$Petal.Length) y <- as_data(iris$Sepal.Length) Species <-as.numeric(iris$Species) n_species <- length(unique(Species)) # variables and priors mu <- normal(0, 10, dim = n_species) sd <- uniform(0, 10, dim = n_species) coef <- normal(0, 10, dim = n_species) # operations mean <- mu[Species] + coef[Species] * x sd_species <- sd[Species] # likelihood distribution(y) <- normal(mean, sd_species) # defining the model m <- model(mu, sd, coef) # plotting plot(m) # sampling draws <- mcmc(m, n_samples = 1000) summary(draws) #result plot mcmc_trace(draws) mcmc_intervals(draws)
今回使ったコード
今回使った全てのコードは以下においてありますので参考にしていただければと思います。
gretaパッケージを活用したベイズモデリングの紹介 · GitHub
所感
かなり簡単にベイズモデリングを行うことが可能なパッケージだと思います。
ただし、一部個人的に物足りなさを感じました。
階層ベイズモデルの書き方がいまいちドキュメントになく書き方が分かりづらい
モデルの可視化をした際に階層性が表現されていない
といった点です。
とはいえ、個人的にはStanよりもかなり簡単にかける気がして、今後積極的に使っていきたいなと感じるほどです。
間違い等ありましたらご指摘お願いいたします。
