Edited at

機械学習でベイジアンモデルを使うという選択肢

More than 1 year has passed since last update.


概要

機械学習やデータ解析を勉強していると時々目にするベイジアンモデルをざっくり理解するために以下の内容をまとめました。


  • どのような背景でベイジアンモデルが研究されているのか

  • イメージをつかむために簡単なベイジアンモデルを実装してみる

  • 一般的な機械学習とベイジアンモデルとの対応


ベイジアンモデルが研究される背景


Zoubin Ghahramani氏の見解


ベイジアンモデルのメリット


  • フレームワークが単純で、一意的。ニューラルネットワークのように任意性のある正則化項が登場したりしない。


  • Gaussian Processのようなノンパラメトリック・ベイジアンモデルは表現力が高い(関数空間のカバー領域が広い)。

  • 過適合が生じない。


ベイジアンモデルのデメリット


  • 並列処理に向かない。

  • モデルの解釈が難しいことがある。

  • アルゴリズムの再現性がない場合がある。

(追記 2017/02/17)

本節の内容は不正確ではないかとのご指摘をいただいています。

本節はこの講義の内容を私なりの解釈でまとめたつもりでしたが、ベイジアンモデル全般に関する主張としては不正確であると思われるので、いただいたご指摘と合わせて元の動画を見ていただければと思います。


機械学習に対する2つの考え方


  • データを色々な手法に入力してみて、良い性能を示すものを使って予測や意思決定をする。

  • モデルが取りうる空間を定義し、データからパラメータを学習して予測や意思決定をする。

どちらのアプローチも一長一短がありますが、ベイジアンモデルは後者の立場に立つ人に好まれます。特に定式化が「エレガントである」点で後者と親和性が高い点が強調されます。


統計モデリング手法として

久保拓弥氏の著書では、現実的な統計モデル構築のためにベイジアンモデルが有用であることが説明されています。例えば空間構造がある階層ベイジアンモデルがわかりやすい応用例です。


簡単なベイジアンモデルの実装


まずは非ベイジアンの統計モデルから

説明変数$x$によって目的変数$y$を説明するモデルを考えます。

このような散布図(=学習データ)から、

non_bayesian_input.png

このように未知の値を予測できる線を引くために、

non_bayesian_regression.png

以下のようなGLMを考えます。

y \sim \mathrm{Pois} (\lambda)

\lambda = \beta_1 + \beta_2 \ln x


入力データ生成

現実問題のデータ解析では以下のように既知の分布から人工的にデータが生成されることはないですが、解析の妥当性を示すために既知の分布を使います。

num_data <- 10

X <- runif(n = num_data, min = 0, max = 2)
beta_1 <- 1
beta_2 <- 3
lambda <- exp(beta_1 + beta_2 * log(X))
Y <- rpois(n = length(lambda), lambda = lambda)


GLMによる回帰

df <- data.frame(x = log(X), y = Y)

glm_result <- glm(y ~ x, data = df, family = "poisson")
df$glm <- glm_result$fitted.values
f <- ggplot(df, aes(x, y))
f + geom_point() + geom_line(aes(x, glm))
glm_result


result

Call:  glm(formula = y ~ x, family = "poisson", data = df)

Coefficients:
(Intercept) x
1.044 3.098

Degrees of Freedom: 9 Total (i.e. Null); 8 Residual
Null Deviance: 72.94
Residual Deviance: 5.081 AIC: 31.59


ポワソン分布のパラメータ$\lambda$を与える切片(1.044)と傾き(3.098)がデータ生成に用いた値$\beta_1 = 1, \beta_2 = 3$にそれなりに近くなりました。


ベイジアンモデルを適用する

このGLMをベイジアンモデルにすることはざっくり言うと「$\beta_1, \beta_2$を確率変数として扱うこと」に対応します。

ベイジアンモデルの大まかな流れは、ベイズの定理

P(パラメータ|目的変数) \propto P(目的変数\vert パラメータ)P(パラメータ)

を用い、


  1. $P(目的変数 \vert パラメータ)$ = 尤度関数を定義(モデリング)し、

  2. $P(パラメータ)$ = 事前分布を定義し、

  3. 何らかの方法(MCMCなど)で$P(パラメータ \vert 目的変数)$ = 事後分布を求め、

  4. 得られた事後分布を予測や意思決定に使う

といったものです。

(この記事は簡便化のために説明変数は定数として扱います)

今回のモデルではデータ数を$n$、$X = \lbrace x_1, \cdots, x_n \rbrace$、$Y = \lbrace y_1, \cdots, y_n\rbrace$とすれば、

P(y_i \vert \beta_1, \beta_2, x_i)

= \frac{\lambda^{y_i}}{y_i! \exp \lambda}
= \frac{(\beta_1 + \beta_2x_i)^{y_i}}{y_i! \exp(\beta_1 + \beta_2 x_i)}

P(Y \vert \beta_1, \beta_2, X) = \prod_i P(y_i \vert \beta_1, \beta_2, x_i)

という関係を利用し、

each_likelihood_function <- function(y, x, beta_1, beta_2) {

dpois(x = y, lambda = exp(beta_1 + beta_2 * log(x)), log = TRUE)
}

likelihood_function <- function(Y, X, beta_1, beta_2) {
each_likelihood_function(Y, X, beta_1, beta_2) %>% sum %>% exp
}

以上の尤度関数を定義しておきます。

$\prod \cdot$の代わりに$\exp (\sum \ln \cdot )$を使うのはアンダーフローを避けるためです。

実装が簡単になるので、以下の$\beta_1, \beta_2$の格子点についての尤度を取得しておきます。

beta_1_candidates <- seq(from = 0.5, to = 1.5, by = 0.1)

beta_2_candidates <- seq(from = 2.5, to = 3.5, by = 0.1)
params_candidates <- merge(x = beta_1_candidates, y = beta_2_candidates, all = TRUE)
names(params_candidates) <- c("beta_1", "beta_2")
likelihood_df <- ddply(params_candidates, .(beta_1, beta_2), function(params) {
data.frame(likelihood = likelihood_function(Y, X, params$beta_1, params$beta_2))
})
filled.contour(
beta_1_candidates,
beta_2_candidates,
likelihood_df$likelihood %>% matrix(nrow = length(beta_1_candidates)),
key.title = title(main = "likelihood"),
xlab = "beta_1", ylab = "beta_2"
)

likelihood_contour.png

GLMや最尤法は、この尤度を最大化させる$\beta_1, \beta_2$を求めていることになります。


無情報事前分布

事前分布$P(\beta_1)$は知ることができず、定義する必要があります。

よく使われるのは「大きな範囲の一様分布」もしくは「ひらべったい正規分布」ですが[1]、この記事では前者を適用して項を無視します。


MCMCGibbs Sampling

尤度関数、事前分布を定義したら、いよいよ事後分布を計算できます。

今回はMCMC(Markov chain Monte Carlo)に分類されるGibbs Samplingを実装してみます。ニューラルネットワークまわりを勉強されている方は、RBMでお馴染みのアルゴリズムです。

$\beta_1, \beta_2$を、一方のパラメータをサンプルするときは他方のパラメータを固定してサンプリングします。事前分布は一様分布であるために無視することができるので、以下の確率分布を用います。

P(\beta_1 \vert Y, \beta_2) \propto P(Y \vert \beta_1, \beta_2)

P(\beta_2 \vert Y, \beta_1) \propto P(Y \vert \beta_1, \beta_2)

この確率はFull conditional distributionといいます。

update_beta <- function(candidates, full_conditional_distribution) {

distribution_function <- function(beta) {
approxfun(x = full_conditional_distribution$beta, y = full_conditional_distribution$likelihood)(beta)
}
threshold <- full_conditional_distribution$likelihood %>% max * 1.1

sample_from_distribution(
candidates = candidates,
distribution = distribution_function,
threshold = threshold
)
}

gibbs_samples <- data.frame(
beta_1 = vector(mode = "numeric", length = num_samples),
beta_2 = vector(mode = "numeric", length = num_samples)
)
beta_1_sample <- head(beta_1_candidates, 1)
beta_2_sample <- head(beta_2_candidates, 1)

for (sample_index in 1:num_samples) {
# update beta_1_sample
full_conditional_distribution_1 <- likelihood_df[likelihood_df$beta_2 == beta_2_sample, c("beta_1", "likelihood")]
names(full_conditional_distribution_1) <- c("beta", "likelihood")
beta_1_sample <- update_beta(
candidates = beta_1_candidates,
full_conditional_distribution = full_conditional_distribution_1
)

# update beta_2_sample
full_conditional_distribution_2 <- likelihood_df[likelihood_df$beta_1 == beta_1_sample, c("beta_2", "likelihood")]
names(full_conditional_distribution_2) <- c("beta", "likelihood")

beta_2_sample <- update_beta(
candidates = beta_2_candidates,
full_conditional_distribution = full_conditional_distribution_2
)
gibbs_samples[sample_index, "beta_1"] <- beta_1_sample
gibbs_samples[sample_index, "beta_2"] <- beta_2_sample
}

ここで、与えられた確率密度関数からサンプリングする手法としてRejection Samplingを使用したので次節に書きます。


Rejection Sampling

詳細はWikipediaをご参照ください。コードだけ載せておきます。

# sample_from_distribution(candidates, distribution ,threshold, n = 1)

# returns a vector of samples that follow a given distribution by the rejection method
# arguments:
# candidates : a vector from which to choose domain samples
# distribution : a distribution function
# threshold : the upper limit for the rejection method
# n : number of samples to return
# usage:
# candidates <- seq(from = -1.0, to = 1.0, by = 0.1)
# distribution <- function(x) dnorm(x)
# threshold <- dnorm(0) * 1.2
# hist(sample_from_distribution(candidates, distribution, threshold))

sample_from_distribution <- function(candidates, distribution, threshold, n = 1) {
sapply(seq(n), function(n_) {
while (TRUE) {
x <- sample(candidates, size = 1)
dice <- runif(n = 1, min = 0, max = 1)
if (dice < distribution(x) / threshold) {
return(x)
}
}
})
}


実行結果

hist(gibbs_samples$beta_1, probability = TRUE)

hist(gibbs_samples$beta_2, probability = TRUE)

を見てみると、以下のヒストグラムを得ます。

beta_1_density.png

beta_2_density.png

生成に用いた真値$\beta_1 = 1, \beta_2 = 3$の周辺の分布を得ることができました。


一般的な機械学習との対応

このようにベイジアンモデルでも結局やることは、与えられたデータセットを学習するために最適化やサンプリングなどの計算を行うという点で他の機械学習と同様です。

概念的には以下ような対応をあてはめるとイメージしやすいです。

一般的な機械学習
ベイジアンモデル

学習
MCMC等

学習パラメータ
確率モデルのパラメータ

評価関数
事後確率


参考文献