R でベイズ線形回帰の予測分布 にて、初 R しつつ、線形回帰の分布図とか書いてみてたら、idojun さんに「バイアス項は？」とつっこまれた。
うああ、忘れてた！
でもバイアス項ってどれくらい効いてくるの？ 無くてもそれっぽい分布得られちゃったんだけど、やっぱり無いとダメなの？ という話になり。
ちょうど PRML 3.4〜5 でモデルの選択を評価するモデルエビデンスを学んだから、いろいろなモデルについてこれを計算してみた。
前回同様、突っ込んでもらいやすいよう step by step で。
前回記事では、基底関数を１パラメータで記述できていたが、そのままガウス＋バイアスに拡張するにはトリッキーな技を駆使する羽目になる。
ここは一般の関数列を扱えるようにしたい。
c() に関数を放り込めば、関数オブジェクトの入ったリストが得られる。
> funcs <- c(function(x)x, function(x)x^2, function(x)x^3) > funcs [[1]] function(x)x [[2]] function(x)x^2 [[3]] function(x)x^3
値を入れて行列を得るには sapply を使う。
> sapply(funcs, function(f)f(1:3)) [,1] [,2] [,3] [1,] 1 1 1 [2,] 2 4 8 [3,] 3 9 27
outer でも、なんかできそうな気がするが、できない。
> outer(funcs, 1:3, function(f,x)f(x)) [,1] [,2] [,3] [1,] 2 3 4 [2,] 2 3 4 [3,] 2 3 4
これで、前回の phi の代わりに任意の関数列を扱う準備ができた。
続いて、いろんな基底関数、PRML に載っている多項式、ガウス、ロジスティック・シグモイドのそれぞれについて「基底関数列を生成する関数」。
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)) } )) }
第１引数は基底関数の次元。
たとえば polynomial_basis_func(3) で 1, x, x^2 の３次元の多項式基底関数が得られる(PRML 本では M=3 は「３次元の多項式」の意味なので注意)。
gaussian_basis_func などの第２引数 has_bias はバイアス項の有無を TRUE/FALSE で指定。
第３引数の s は「空間的な尺度」とやら。
PRML本には s の決め方について一言も書いていないので、適当でいいのかと思っていたら、漠然と危惧していたよりもずっと敏感に効いてくるようだ。が、それはひとまずおいておく。
なお。
バイアス項を記述するのに "function(x) 0*x+1" となっていて、「なんで "function(x) 1" じゃあダメなの？」とかいうあたりは、「 R ハマリどころ」とか別のまとめ記事とか書こうかな。
実はさっきの「 outer でもできそうなのに、できない」のも同じ理由。
こうして準備した基底関数列を元に、モデルのエビデンスを計算する。
p(M|D) = p(α,β|D) ∝ p(D|α,β)p(α,β)
PRML 3.5.1 では、ln p(D|α,β) を求めて、これが「エビデンスの表式」だと言うんだけど、モデルの事前分布 p(α,β) はどこいった？
事前分布は均等ということにでもして、定数比例の範囲で無視してんのか？
強いて偏りがあるとすればβの真値だろうけど、それを知ってるなんて前提嬉しくないし。
まあでも「エビデンスの表式」なんだから ln p(D|α,β) を計算すればよいということになっているのだろう。
PRML 本の (3.81)(3.84)(3.82)、そして (3.86) から ln p(D|α,β) を求められる。
# alpha, beta, D, phi はあらかじめ与える M <- length(phi) N <- length(D$x) PHI <- sapply(phi, function(f)f(D$x)) # design matrix A <- alpha * diag(M) + beta * t(PHI) %*% PHI # (3.81) equal to S_N (3.54) m_N <- beta * solve(A) %*% t(PHI) %*% D$t # (3.84) loss_m_N <- D$t - PHI %*% m_N E_m_N <- beta / 2 * sum(loss_m_N^2) + alpha / 2 * sum(m_N^2) # (3.82) evidence <- M/2*log(alpha) + N/2*log(beta) - E_m_N - 1/2*log(det(A)) - N/2*log(2*pi) # (3.86)
数式をそのまま記述できるのはＲの嬉しいところ（いろいろコツはいるけど）。
まとめて、利用しやすいように関数の形にして、機能をちょいちょい足すとこんな感じ。
βの最尤値 求めるとか、グラフを描くとか、機能を増やしていくうちにえらく長くなってしまった……
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) }
ハイパーパラメータのαは明示的に与える必要があるが、βは "ml" を指定すれば最尤値(PRML 1.63式 や 3.21式)を求めてそれを用いる。
本当は PRML 3.5.2 の方法でαとβを求めるべきなのだろうけど、ちょっとめんどくさそうなので、とりあえずパス。
さて。
こうして実装したエビデンス計算が正しいかどうか確認したい。
PRML p167 の図 3.14 は、多項式基底関数の次数に対してモデルエビデンス(の対数)がどのような値を取るかを表している。
データセットはサポートサイトで公開されており、ハイパーパラメータαも本文中に 5×10^-3 と指定があるので、それらを使って多項式基底関数のエビデンスを計算すれば、このグラフ通りになることが期待される。
唯一、βだけが不明だが、最尤値か真の値 1/0.3^2 = 11.1 のどちらかではないか、と推察。
というわけで２種類のβに対してエビデンスを求めたのが以下の結果。
次元 M はグラフと見比べられるよう、基底関数の個数ではなく多項式の次数とした(なので基底関数の個数は M+1)
# PRML's synthetic data set curve_fitting <- data.frame( x=c(0.000000,0.111111,0.222222,0.333333,0.444444,0.555556,0.666667,0.777778,0.888889,1.000000), t=c(0.349486,0.830839,1.007332,0.971507,0.133066,0.166823,-0.848307,-0.445686,-0.563567,0.261502)) > # ０次〜８次の多項式基底関数のエビデンスを求める > # βは最尤推定値の場合 > a<-sapply(1:9, function(n) calc_evidence(polynomial_basis_func(n), curve_fitting, alpha=5e-3, beta="ml")) > data.frame(M=0:8, evidence=a[1,], beta_ML=a[2,]) M evidence beta_ML 0 -13.60463 2.649926 1 -14.48098 4.680463 2 -16.60761 4.752649 3 -14.38654 28.600038 4 -14.20562 28.651286 5 -15.12706 29.206330 6 -15.86874 30.294868 7 -16.43925 30.954700 8 -17.37590 35.353486 > # 真の値 β = 1/0.3^2 = 11.1 を用いる場合 > a<-sapply(1:9, function(n) calc_evidence(polynomial_basis_func(n), curve_fitting, alpha=5e-3, beta=11.1)) > data.frame(M=0:8, evidence=a[1,], beta=a[2,]) M evidence beta 0 -23.10268 11.1 1 -17.88419 11.1 2 -20.30879 11.1 3 -13.93411 11.1 4 -13.71294 11.1 5 -14.35868 11.1 6 -14.98120 11.1 7 -15.48112 11.1 8 -15.90587 11.1
ん〜〜〜……ビミョウ。
β＝11.1 のときは結構いい線をいってる。ただし、僅差で M=3 より M=4 のエビデンスの方が大きい。
最尤値の場合は M=0 がエビデンス最大（ええええ？）。
なんでこんなことになるのか、式を見ながら考える。
M=0 だと点を外しまくっているわけだから、当然精度パラメータβの最尤値は小さくなる。M=3 だといい感じの回帰関数が得られるので精度の最尤値は大きくなる。
一方 (3.82) から、精度βが高い場合は、点を外してしまった場合の損失関数への影響が大きくなる。
実際計算してみると、M=0 より M=3 の時の方が E(m_N) は大きくなってしまっており、これがエビデンスの大小にそのまま効いてきていた。
これは、βを動かしながらモデル間のエビデンスを比較することはできないということを指しているのかなあ。うーん。
ちょっと長くなっちゃったので、本題だった「バイアス項の有無でエビデンスはどう変わるか？」は次回乞うご期待で。
→続き 「基底関数を色々変えて、線形回帰のエビデンスを計算してみた part 2」
