概要

『ベイズ推論による機械学習入門』のChapter3のp105~p114に線形回帰の事後分布の求め方がのっていたので、実際にJuliaで実装し、線形回帰のパラメータを求めてみました。

機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)

作者: 須山敦志,杉山将
出版社/メーカー: 講談社
発売日: 2017/10/21
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

ちなみに今回のアプローチは経験ベイズ的なアプローチです。

完全ベイズ的なアプローチである、事前分布のハイパーパラメータに自体に分布を仮定してMCMCでサンプリングして更新していくというやり方ではなく、事前分布のパラメータは適当な値を設定し*1、データから計算して事後分布のパラメータを求めていくアプローチです*2。

したがって既にStanなどで完全ベイズ的なベイズ推定に親しんでいて、そちらしか知らない方には少々違和感を感じる部分があるかもしれませんが、本稿では経験ベイズ的なアプローチとなり異なりますのでご注意ください。

今回使った全体のコードは以下です。

bayes_linear_regression.ipynb · GitHub

モデルの構築

実装に入る前に本書での数式を記述しておきます。

線形回帰モデル

線形回帰モデルは以下のようにモデル化されます。

${y_{n} = \textbf{w}^{T} \textbf{x}_{n} + e_n}$

ここで、y_nは目的変数、x_nが説明変数、wは重みパラメータe_nがノイズ（誤差）です。

ノイズ成分は以下のように平均ゼロのガウス分布と仮定。

${e_{n} \sim N(e_n |0, \lambda^{-1})}$

ここで出てきたλは既知の精度パラメータとします*3。

確率分布として定式化すると、

${p(y_n|\textbf{x}_{n}, \textbf{w}) = N(y_n | \textbf{w}^{T} \textbf{x}_{n}, \lambda^{-1}) }$

となります。

事前分布と事後分布に関して

今回はパラメータwを観測したデータから学習するために事前分布は以下のような正規分布を仮定します。

${p(\textbf{w}) = N(\textbf{w} | \textbf{m}, \Lambda^{-1})}$

このmは平均パラメータで、Λは正定値行列の精度行列パラメータとなります。これらはハイパーパラメータとなり事前に決める必要があります。

一方でwの事後分布の方はというとベイズの定理を利用して、

${ p(\textbf{w}| \textbf{Y}, \textbf{X}) = \frac{p(\textbf{w}) \Pi^{N}_{n = 1} \ p(y_{n} | \textbf{x}_{n}, \textbf{w} )} {p(\textbf{Y}| \textbf{X})}\\ \propto p(\textbf{w}) \Pi^{N}_{n = 1} p(y_{n} | \textbf{x}_{n}, \textbf{w} ) }$

と表現できます。

この式に対して対数に入れてwで整理するとパラメータ計算のためのが求まります*4。

最終的には事後分布は以下の式になります。

${ p(\textbf{w}| \textbf{Y}, \textbf{X}) = N(\textbf{w} | \hat{\textbf{m}}, \hat{\Lambda}^{-1}) }$

ただし、

${ \hat{\Lambda}^{-1} = \lambda \sum^{N}_{n=1} \textbf{x}_{n} \textbf{x}_{n}^{T} + \Lambda\\ \hat{\textbf{m}} = \hat{\Lambda}^{-1}(\lambda \sum^{N}_{n=1} y_n \textbf{x}_{n} + \Lambda \textbf{m}) }$

となり、上のm_hatとΛ_hatはデータから計算し、その計算したパラメータからwを生成すれば良いことがわかります。

事後分布のパラメータを求める関数の実装としては以下のようにしています。muとlambdaは特に指定しない場合は勝手に値が入るようにしました*5。

function weight_distribution_hyper_parameter(y, X, mu, lambda, l)
    if mu == 0 & lambda == 0
        #初期値
        mu = zeros(size(X)[1])
        lambda = eye(size(X)[1]) #単位行列
    end

    lambda_hat = (l * X * X') + lambda
    lambda_hat_inv = inv(lambda_hat)

    mu_hat = lambda_hat_inv * ((l * y' * X') + mu' * lambda)'
    
    return [mu_hat, lambda_hat, lambda_hat_inv]
end

この関数は実際には自分で以上の関数を直接使うことは想定していなく、以下の学習の関数の中で呼ぶようにしています。

4行目でパラメータを計算して、6行目でwをサンプリングしています。

インターフェースは若干scikit-learn風で、各種パラメータを返します。

function bayes_fit(formula, df; mu = 0, lambda = 0, l = 200, iter = 100)
    
    dm = ModelMatrix(ModelFrame(formula, df)).m' #計画行列
    parameter = weight_distribution_hyper_parameter(df[formula.lhs], dm, mu, lambda, l) #パラメータ計算結果
    
    weight = sampling_weight(parameter[1], parameter[3], iter)
    
    return [parameter, weight, formula, df, dm, l]
end

予測分布に関して

予測分布に関しては詳細は省略しますが、以下の式で表現できます*6。

${ p(y_{*}| \textbf{x}_{*}) = N(y_{*}| \mu_{*}, \lambda^{-1}_{*} )\\ }$

ただし、

${ \mu_{*} = \textbf{m}^{T} \textbf{x}_{*} \\ \lambda^{-1}_{*} = \lambda^{-1} + \textbf{x}_{*}^{T} \Lambda^{-1} \textbf{x}_{*} }$

となります。このmとΛをデータから求めた事後分布のパラメータであるm_hatとΛ_hatで置き換えれば、データから学習した結果で予測分布を生成できます。

実装上は以下のようにしました。

function predict(model, X)
    df = deepcopy(X)
    df[model[3].lhs] = 1 #目的変数に適当な値を設定、特に意味はない
    df = ModelMatrix(ModelFrame(model[3], df)).m'
    
    mu_p = (model[1][1]' * df)'
    l_p = diag(1 / model[6] + df' * model[1][3] * df) #対角が精度
    
    mu_p_plus = mu_p + sqrt.(l_p)
    mu_p_minus = mu_p - sqrt.(l_p)
    
    return [mu_p mu_p_plus mu_p_minus l_p]
end

こちらも学習の関数と同様にscikit-learn風で、学習関数の返り値を第一引数に、第二引数に予測したいデータを入れます。

返り値は予測値と予測値に±√λ^-1した値です。

モデルエビデンスに関して

モデルエビデンスとはベイズ学習において複数のモデルを比較する際に使われる値です。

計算式としては事後分布の分母（p(Y|X)）の値の対数をとったものになります*7。したがって、

${ ln \ p(\textbf{Y}| \textbf{X}) = -\frac{1}{2} \{ \sum^{N}_{n=1}(\lambda y^{2}_{n} - ln \ \lambda + ln \ 2 \pi) + \textbf{m}^{T} \Lambda \textbf{m} - ln \ |\Lambda| - \hat{\textbf{m}}^{T} \hat{\Lambda} \hat{\textbf{m}} + ln \ |\hat{\Lambda}| \} }$

で求めることができます。

実装は以下のようにしました。

#モデルエビデンス計算
function model_evidence(model)
    #各種パラメータ
    mu = zeros(size(model[5])[1])
    lambda = eye(size(model[5])[1]) #単位行列
    mu_hat = model[1][1]
    lambda_hat = model[1][2]
    y = model[4][:y]
    l = model[6]  

    me = -0.5 * (sum((l * (y .^ 2)) - log(l) + log(2 * pi)) + mu' * lambda * mu - log(det(lambda)) - mu_hat' *  lambda_hat * mu_hat + log(det( lambda_hat)))
    
    return me
end

こちらも引数は予測分布の関数と似ていて、学習関数の返り値を入れるようにしました。

実験

最後にこれまで見てきた事後分布、予測分布の確率計算式の結果を活用して適当なデータで学習してみたいと思います。

流れとしては適当なデータを生成し複数パターンの次数の多項式モデルを学習して、モデルエビデンスの値を見ながら比較をするというものになります。

まずはデータの生成です。

#seed
srand(1234)

#データ生成
x = rand(Uniform(-1, 7), 10)
y = sin.(x)
df = DataFrame(x = x, y = y)

#可視化
plot(df, x=:x, y =:y, Geom.point)

こんなデータが得られました。

f:id:rf00:20180809213714p:plain

フィットするためのデータセットの準備、フィット、学習データと同じ関数で生成した新規データの予測、そして予測結果の可視化をしていきます。

#各種モデル用データ準備
df2 = DataFrame(x = x, x2 = x.^2, y = y)
df3 = DataFrame(x = x, x2 = x.^2, x3 = x.^3, y = y)
df4 = DataFrame(x = x, x2 = x.^2, x3 = x.^3, x4 = x.^4, y = y)
df5 = DataFrame(x = x, x2 = x.^2, x3 = x.^3, x4 = x.^4, x5 = x.^5, y = y)
df10 = DataFrame(x = x, x2 = x.^2, x3 = x.^3, x4 = x.^4, x5 = x.^5, x6 = x.^6, x7 = x.^7, x8 = x.^8, x9 = x.^9, y = y)

#モデルフィット
model0 = bayes_fit(@formula(y ~ 1), df, l = 200, iter = 1000)
model1 = bayes_fit(@formula(y ~ x), df, l = 200, iter = 1000)
model2 = bayes_fit(@formula(y ~ x + x2), df2, l = 200, iter = 1000)
model3 = bayes_fit(@formula(y ~ x + x2 + x3), df3, l = 200, iter = 1000)
model4 = bayes_fit(@formula(y ~ x + x2 + x3 + x4), df4, l = 200, iter = 1000)
model5 = bayes_fit(@formula(y ~ x + x2 + x3 + x4 + x5), df5, l = 200, iter = 1000)
model10 = bayes_fit(@formula(y ~ x + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9), df10, l = 200, iter = 1000)

#seed
srand(2234)

#新規データ予測
x_new = rand(Uniform(-1, 6), 100)

df_new0 = DataFrame(x = x_new)
df_new2 = DataFrame(x = x_new, x2 = x_new.^2)
df_new3 = DataFrame(x = x_new, x2 = x_new.^2, x3 = x_new.^3)
df_new4 = DataFrame(x = x_new, x2 = x_new.^2, x3 = x_new.^3, x4 = x_new.^4)
df_new5 = DataFrame(x = x_new, x2 = x_new.^2, x3 = x_new.^3, x4 = x_new.^4, x5 = x_new.^5)
df_new10 = DataFrame(x = x_new, x2 = x_new.^2, x3 = x_new.^3, x4 = x_new.^4, x5 = x_new.^5, x6 = x_new.^6, x7 = x_new.^7, x8 = x_new.^8, x9 = x_new.^9)

#新規データ予測
pre0 = predict(model0, df_new0)
pre1 = predict(model1, df_new0)
pre2 = predict(model2, df_new2)
pre3 = predict(model3, df_new3)
pre4 = predict(model4, df_new4)
pre5 = predict(model5, df_new5)
pre10 = predict(model10, df_new10)

#予測と実績を可視化
set_default_plot_size(30cm, 30cm)#図のサイズ設定

p0 = plot(layer(x=x, y=y, Geom.point, Theme(default_color="black")),
layer(x=x_new, y=pre0[:, 1], Geom.line, Theme(line_style =:solid, default_color="black")),
    layer(x=x_new, y=pre0[:, 2], Geom.line, Theme(line_style =:dash)),
    layer(x=x_new, y=pre0[:, 3], Geom.line), Theme(line_style =:dash))

p1 = plot(layer(x=x, y=y, Geom.point, Theme(default_color="black")),
layer(x=x_new, y=pre1[:, 1], Geom.line, Theme(line_style =:solid, default_color="black")),
    layer(x=x_new, y=pre1[:, 2], Geom.line, Theme(line_style =:dash)),
    layer(x=x_new, y=pre1[:, 3], Geom.line), Theme(line_style =:dash))

p2 = plot(layer(x=x, y=y, Geom.point, Theme(default_color="black")),
layer(x=x_new, y=pre2[:, 1], Geom.line, Theme(line_style =:solid, default_color="black")),
    layer(x=x_new, y=pre2[:, 2], Geom.line, Theme(line_style =:dash)),
    layer(x=x_new, y=pre2[:, 3], Geom.line), Theme(line_style =:dash))


p3 = plot(layer(x=x, y=y, Geom.point, Theme(default_color="black")),
layer(x=x_new, y=pre3[:, 1], Geom.line, Theme(line_style =:solid, default_color="black")),
    layer(x=x_new, y=pre3[:, 2], Geom.line, Theme(line_style =:dash)),
    layer(x=x_new, y=pre3[:, 3], Geom.line), Theme(line_style =:dash))


p4 = plot(layer(x=x, y=y, Geom.point, Theme(default_color="black")),
layer(x=x_new, y=pre4[:, 1], Geom.line, Theme(line_style =:solid, default_color="black")),
    layer(x=x_new, y=pre4[:, 2], Geom.line, Theme(line_style =:dash)),
    layer(x=x_new, y=pre4[:, 3], Geom.line), Theme(line_style =:dash))

p5 = plot(layer(x=x, y=y, Geom.point, Theme(default_color="black")),
layer(x=x_new, y=pre10[:, 1], Geom.line, Theme(line_style =:solid, default_color="black")),
    layer(x=x_new, y=pre10[:, 2], Geom.line, Theme(line_style =:dash)),
    layer(x=x_new, y=pre10[:, 3], Geom.line), Theme(line_style =:dash))


gridstack([p0 p1 p2; 
    p3 p4 p5])

各データセット毎に以下のような予測結果と予測の幅が得られました。

f:id:rf00:20180809214816p:plain

左上から右下に行くに連れて複雑なモデルで予測していて、だんだんと精度がよくなっている感じが見えています。

ただし一番複雑な右下のモデルは幅が広くなっている部分があり、不確実さが増している感じがあります。

モデルエビデンス観点から比較してみると、

#モデルエビデンス可視化
me = zeros(2, 6)
me[1, 1] = model_evidence(model0)
me[1, 2] = model_evidence(model1)
me[1, 3] = model_evidence(model2)
me[1, 4] = model_evidence(model3)
me[1, 5] = model_evidence(model4)
me[1, 6] = model_evidence(model10)

me[2, 1] = 1
me[2, 2] = 2
me[2, 3] = 3
me[2, 4] = 4
me[2, 5] = 5
me[2, 6] = 10

set_default_plot_size(10cm, 10cm)#図のサイズ設定
plot(x=me[2, :], y=me[1, :], Geom.line)

4つ目の左下のモデルが一番良さそうです。

左下と真ん中下は一見どちらも良さそうですが、詳しく見てみるとxが小さいときの幅が真ん中で大きくなっていて、多少不確実さが高くなっている様子が見られました。

f:id:rf00:20180809215405p:plain

まとめ

今回事前分布から事後分布を解析的に求めた式を追いながらコードに落とし込むということをやってみてとても勉強になりました。

Stanでベイズ推定をやっていただけでは身につかないベイズの奥深さを少し理解できた気がします。

間違いなどありましたら指摘お願いいたします。

*1:事前知識を盛り込むこともできる

*2:詳しくは『ベイズ推論による機械学習入門』のp79などを参照のこと

*3:実際に現実世界で既知であることはあるのか個人的にはよくわかってないです。コードの実装上は任意に好きな値を入れられるようにしているのですが、どういうふうに決めたらいいのかご存知の方いたら教えてください

*4:ここでは略しますが、『ベイズ推論による機械学習入門』のp107の式(3.146)です

*5:何らかの事前知識がある場合は入れてもよいはず

*6:『ベイズ推論による機械学習入門』のp107, p108に説明があります

*7:こちらも詳細は『ベイズ推論による機械学習入門』のp109, p110に説明があります

Deep Karmaning

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

『ベイズ推論による機械学習入門』の線形回帰の例をJuliaで実践

概要

モデルの構築

線形回帰モデル

事前分布と事後分布に関して

予測分布に関して

モデルエビデンスに関して

実験

まとめ