Deep Karmaning

技術系の話から日常のことまで色々と書きます

gretaパッケージを活用したベイズモデリングの紹介

概要

今回はRのgretaパッケージを使ったベイズモデリングを紹介します。

公式サイトは以下。

https://greta-dev.github.io/greta/index.html

公式のサンプルコードは以下にのっています。

Get started with greta

gretaとは

gretaパッケージはどんなパッケージかと言うと、tendorflowをバックで動かしてモデリング可能なパッケージです。

色々出来るようなのですが、今回は公式サイトでも紹介されているようなMCMCを用いたベイズモデリングを試します。

特徴としては以下の3点があげられていて、

  • simple

    • Rで書かれていて、BUGSやStanを覚える必要がない
  • scalable

    • Google TensorFlowを使っているので大規模なデータにも対応可能、CPUでもGPUでも動く
  • extensive

    • R関数を簡単にかける、また他のRパッケージを活用できる

ということらしいです。

これだけ読むと確率的プログラミング言語として、BUGSやStanが意識されている気がします。

実際、サイトのexampleを読むと、BUGSとStanのコードと比較した物が乗っており、gretaでは短くシンプルに書くことが出来ることを示されています。

インストール方法

まずは以下記事を参考にgretaパッケージを入れます。

Get started with greta

install.packages("greta")

次にtensorFlowが必要なので、インストールしましょう。

インストール用の関数はgretaパッケージを読み込むと用意されていますので、以下を実行。

library(greta)

install_tensorflow()

また以下のようなお洒落な図を書くためにDiagrammeRを入れます(図の意味は後述)。

f:id:rf00:20180506130630p:plain

install.packages('DiagrammeR')

これでgretaを使う準備ができました。

モデリング

今回はモデリングではirisデータを使って二つのモデルを試したいと思います。

具体的には、

  1. Sepal.Lengthの平均とばらつきを種ごとに推定するモデル

  2. 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)

こういった図が表示されると思います。

f:id:rf00:20180506145415j:plain

この図の各図形の意味は公式サイトにも乗っていますが、

f:id:rf00:20180506145542p:plain

  • 白抜きの四角:データ

  • 灰色の丸:変数

  • ひし形:確率分布

  • 小さい丸:演算

  • 線:決定した値

  • 点線:確率値

となります。

先ほどプロットした図では種ごとに異なるといった階層性がうまく表現されていないので少し分かりづらいですが、 実際には種ごとに平均とばらつきが異なるといったモデリングしています。

図は階層性が表現されていないことを目をつぶるととてもわかり易いですね。

数式で表現すると、

{y 〜 Normal(\mu[k], sd[k])}

{\mu[k] 〜 Normal(0, 10)}

{sd[k] 〜 Uniform(0, 10)}

なります。

ここで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)

f:id:rf00:20180506152302p:plain

f:id:rf00:20180506152358p:plain

ちなみに、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よりもかなり簡単にかける気がして、今後積極的に使っていきたいなと感じるほどです。

間違い等ありましたらご指摘お願いいたします。