前回長くなりすぎて力尽きた「線形回帰のエビデンス」の続き。って今回も むやみに長いんだけど。
パターン認識と機械学習(PRML)読書会 #4 で idojun さんと「バイアス項って要るの？」という話になり。
PRML 3.4〜3.5 のモデルエビデンス使えば評価できるかも！ きっとバイアス項有りの方がエビデンスが大きくなるんだよ！ ＯＫ、計算してみた！ というお話。
もちろん結果はとっくに出てたんだけど、なかなか記事に仕立てる間が無く。
録りたまってた東のエデンもようやく見終えた*1ので、そろそろ書く。
基底関数
まずは多項式、ガウス、ロジスティック・シグモイドの各基底関数列を生成する関数(再掲)。
第１引数は基底関数の次元。
第２引数 has_bias はバイアス項の有無を TRUE/FALSE で指定。
第３引数の s は「空間的な尺度」とやら。
polynomial_basis_func <- function(M) { lapply(1:M, function(u){ u1=u; function(x) x^(u1-1) }) } gaussian_basis_func <- function(M, has_bias=T, s=0.1) { phi <- c() if (has_bias) phi <- function(x) 0*x+1 # bias u_i <- seq(0,1,length=ifelse(has_bias, M-1, M)) append(phi, lapply(u_i, function(u){ u1=u; function(x) exp(-(x-u1)^2/(2*s*s)) })) } sigmoid_basis_func <- function(M, has_bias=T, s=0.1) { phi <- c() if (has_bias) phi <- function(x) 0*x+1 # bias u_i <- seq(0,1,length=ifelse(has_bias, M-1, M)) append(phi, lapply(u_i, function(u){ u1=u; function(x) 1/(1+exp(-(x-u1)/s)) } )) }
エビデンス関数
次に、モデルエビデンスを求める関数(再掲)。
第１引数は基底関数列、第２引数は訓練データ。
第３＆４引数はハイパーパラメータのαとβ。αは明示的に与える必要があるが、βは "ml" を指定すれば最尤値(PRML 1.63式 や 3.21式)を求めてそれを用いる。
第５引数に真の関数を与えると、グラフを描く。
calc_evidence <- function(phi, D, alpha=2, beta=25, graph=NULL) { M <- length(phi) N <- length(D$x) PHI <- sapply(phi, function(f)f(D$x)) if (!is.null(graph)) { plot(graph, lty=2, col="blue", xlim=c(0,1), ylim=c(-1.1,1.1), ylab="") par(new=T) plot(D, xlim=c(0,1), ylim=c(-1.1,1.1), xlab="", ylab="") } # maximum likelyhood if (beta=="ml") { w_ML <- solve(t(PHI) %*% PHI) %*% t(PHI) %*% D$t loss_ML <- D$t - PHI %*% w_ML beta <- N / sum(loss_ML^2) if (!is.null(graph)) { par(new=T) plot( function(x) sapply(phi, function(f)f(x)) %*% w_ML , col="red", xlim=c(0,1), ylim=c(-1.1,1.1), ylab="") } } A <- alpha * diag(M) + beta * t(PHI) %*% PHI # equal to S_N(PRML 3.54) m_N <- beta * solve(A) %*% t(PHI) %*% D$t loss_m_N <- D$t - PHI %*% m_N E_m_N <- beta / 2 * sum(loss_m_N^2) + alpha / 2 * sum(m_N^2) if (!is.null(graph)) { par(new=T) plot( function(x) sapply(phi, function(f)f(x)) %*% m_N, xlim=c(0,1), ylim=c(-1.1,1.1), ylab="") } # model evidence c(M/2*log(alpha) + N/2*log(beta) - E_m_N - 1/2*log(det(A)) - N/2*log(2*pi), beta) }
訓練データ。
PRML のサポートサイトで配布されているデータは 10 点しか無くてちょっと役不足なので、お手製の sin 2πx + N(0, 0.2) (25点)を用意。
xlist <- seq(0, 1, length=25) tlist <- sin(2 * pi * xlist)+rnorm(length(xlist), sd=0.2) D <- data.frame(x=xlist, t=tlist)
シグモイド基底関数のモデルエビデンス
９次元のシグモイド関数について、バイアス有無でモデルエビデンスを計算してみる。
結果の１項目がエビデンス、２項目がβ(最尤値)。
負の値なので一瞬惑うが、「バイアス項有りの方がエビデンスが大きい」。
> # with bias > calc_evidence(sigmoid_basis_func(9, T), D, alpha=2, beta="ml", function(x)sin(2 * pi * x) ) [1] -11.37790 41.79888 > # witout bias > calc_evidence(sigmoid_basis_func(9, F), D, alpha=2, beta="ml", function(x)sin(2 * pi * x) ) [1] -12.45082 41.48465
実は、ロジスティック・シグモイド関数の形(しかも位置を表す値 μ_i は、定義域の均等割)から、バイアス項がないと x=0 付近での「解像度」（という言い方は変かな〜？）が落ちてしまうことが推察されるので、バイアス項無しのモデルはもともと適していないはずだったりする。
その意味では、正しく「適したモデルかどうか」を反映した値になっている
ちなみにバイアス有りの場合のグラフ。
青がガウス分布による線形回帰。赤が最尤推定、点線が真のグラフ。
こちらがバイアス無し。他は上とそっくりだが、x=0 付近だけちょっと外してしまっている。
ガウシアン基底関数のモデルエビデンス
ほかのモデルでも比較してみよう。
９次元のガウス基底関数で、αは 2 とし(PRML 3.1)、βは最尤値を用いた。
「バイアス項無しの方がエビデンスが大きい」。
> # with bias > calc_evidence(gaussian_basis_func(9, T), D, alpha=2, beta="ml") [1] -6.62069 41.90476 > # witout bias > calc_evidence(gaussian_basis_func(9, F), D, alpha=2, beta="ml") [1] -6.265361 43.413693
同じ９次元同士だと、「ガウシアン関数８個＋バイアス」と「ガウシアン関数９個」を比べていることになり、バイアス有りの方が不利な気がする。
でも「ガウシアン関数８個＋バイアス」と「ガウシアン関数８個」を比べるのは、今度は逆にバイアス有りの方が断然有利な気がするが、一応計算してみた。
「バイアス項無しの方がエビデンスが大きい」。
> # with bias > calc_evidence(gaussian_basis_func(9, T), D, alpha=2, beta="ml") [1] -6.62069 41.90476 > # without bias > calc_evidence(gaussian_basis_func(8, F), D, alpha=2, beta="ml") [1] -5.890924 41.761091
前回、βに最尤値を用いると、１次の多項式基底関数がエビデンス最大という頭の痛い結果になり、真の値を用いると、そこそこ悪くない結果だった。
そこでβに真の値 1/0.2^2 = 25 を指定してみたのだが、それでも「バイアス項無しの方がエビデンスが大きい」。
> calc_evidence(gaussian_basis_func(9, T), D, alpha=2, beta=25) [1] -6.035784 25.000000 > calc_evidence(gaussian_basis_func(9, F), D, alpha=2, beta=25) [1] -5.585562 25.000000
「空間的な尺度」
実は、バイアスの有無より遙かにエビデンスに効いてくるのが「空間的な尺度」。
gaussian_basis_func() の第３引数を適当に変えると、s=0.2 のあたりに極大がありそうなことがわかる。
> calc_evidence(gaussian_basis_func(9, T, 0.05), D, alpha=2, beta="ml", function(x)sin(2 * pi * x) ) [1] -12.26722 36.55397 > calc_evidence(gaussian_basis_func(9, T, 0.1), D, alpha=2, beta="ml", function(x)sin(2 * pi * x) ) [1] -6.62069 41.90476 > calc_evidence(gaussian_basis_func(9, T, 0.2), D, alpha=2, beta="ml", function(x)sin(2 * pi * x) ) [1] -3.867634 41.629616 > calc_evidence(gaussian_basis_func(9, T, 0.3), D, alpha=2, beta="ml", function(x)sin(2 * pi * x) ) [1] -16.41628 40.93659
s=0.05 のグラフ。
s=0.1 のグラフ。
s=0.2 のグラフ。
エビデンスの大小とモデルの良否が対応している、ように思える。
結論
いろいろなことを推定するのも ためらわれる程度の結果ではあるけど、
- モデルエビデンスは、モデルの評価をある程度表せている
- バイアスの有無が一概にモデルの良否を決めるわけではない
という「雰囲気」は感じられるのでは、と思う。
余談
βの最尤値が真の値 (1/0.2^2 = 25) とかけ離れてるな〜、というのは、点の数を増やせばＯＫ。
> xlist <- seq(0, 1, length=200) > tlist <- sin(2 * pi * xlist)+rnorm(length(xlist), sd=0.2) > D200 <- data.frame(x=xlist, t=tlist) > calc_evidence(gaussian_basis_func(9), D200, alpha=2, beta="ml", function(x)sin(2 * pi * x) ) [1] 16.46175 25.53402
グラフも真の関数に近づいていく。
*1:続きは劇場で、って どゆこと！？
