R でベイズ線形回帰の予測分布

R 機械学習 PRML

一番は「やっぱりＲ覚えよう……」としみじみ実感したことかもしれない(苦笑)。

というわけで R 始めました。


同じことやっても仕方ないので PRML 3.3.2 のベイズ線形回帰による予測分布をやってみることに。
とはいえ、昨日インストールして、今日 R-Tips を読みながら使い始めたという、おしりに立派な殻がついた ひよこちゃんなので、心優しい人が突っ込んでくれることを期待して、 step by step で書いてみる。


R-Tips を読んで、R とは

  • ベクトルと行列の操作が全ての基本。とても便利に使えるようになっている
  • ベクトルとスカラーを演算すると、スカラーが自然にベクトルに拡張される。たとえば v + 3 は ベクトル v の全要素に 3 を加えたベクトルを返す
  • スカラーを取る関数にベクトルを与えることもできる。たとえばベクトル (1,2,3) を log() に入れると、(log(1), log(2), log(3)) を返す。
  • d[分布名] で確率密度関数を、r[分布名] でその分布に従った乱数を返す。主立った確率分布は一通り揃っている(当然)

という感じのものと理解。とりあえずこれでサンプルソースとかが読めるようになった。



ではまずは手始めのウォーミングアップに、１２個の乱数で本当に正規分布の近似になるの？　をヒストグラムを描いて ふむふむしてみる。


乱数を返す関数ということで、rand() とか rnd() とか random() とか探しても見つからない。


r[分布名]() が求めるもの。
いわゆる「単なる乱数」は一様分布(unif)から生成される乱数なので、

> runif(1)
[1] 0.2529997

この乱数が 12個だからって、runif(1)+runif(1)+…… ってエクセルみたいにするんじゃあなくて

> runif(12)
 [1] 0.5778066 0.7346547 0.7584891 0.3959031 0.1103321 0.8387258 0.1636878
 [8] 0.5183256 0.1198709 0.1031732 0.5930027 0.6990498

これが R での「ベクトル」。
横に 12 個並べられないので改行表示。"[8]" は２行目の先頭のインデックス。
すると runif(1) の結果もよく見るとわかるが、実は純粋なスカラーではなく、なんと実は１次のベクトルだったのだ！（←びっくりするところ）


ベクトルの和を取るには sum() が使える。

> sum(runif(12))-6
[1] 0.5424602

６引いて(きっと)正規分布 N(0, 1) に従う(はずの)乱数が１つ得られた。


ヒストグラムを描くには、これを 10000個くらいは用意したい。
この程度で for 文とか使ってしまっては、にわかＲerの名折れ。
しばらく考える。


こうかな？

> rowSums(matrix(runif(12*20), ncol=12))-6
 [1]  1.0596418 -0.6493426 -0.2594199  1.2228629 -0.3831486 -1.1665972
 [7] -0.4137170  0.6124912 -0.3443323 -1.1196722  1.8745801 -0.1115903
[13]  0.2533000  0.6238234 -0.2254362  0.9498943 -1.0495125  1.1735328
[19] -0.6846950  1.5077923

matrix(runif(12 * N), ncol=12) で 12×(欲しい個数) の行列を作って乱数で埋めてから、行ごとに和を取っている。
なんかむやみに富豪的で、ちょっぴり気が咎めるけど。


あとはヒストグラムに突っ込めばＯＫ。

hist(rowSums(matrix(runif(12*10000), ncol=12))-6, breaks=30)

f:id:n_shuyo:20190427122000p:plain

できたできた。きれいにできた。
って、１行で書けちゃうんかいな！



では本題。
ベイズ線形回帰の予測分布。


lm (線形回帰パッケージ)は御法度、の縛りプレイ。
Ｒのお勉強以前に、線形回帰のお勉強なので、まんまのパッケージ使っちゃあねえ。
結論から言うと、使わなくても短く書けるので大丈夫。


線形回帰には基底関数(basis function)が必要。
PRML 3.3.2 ではガウス基底関数 \phi_j(x) = \exp \left{-\frac{(x-\mu_j)^2}{2s^2}\right} を使っているので、まずはこれを記述。

> s <- 0.1
> u_i <- 0.5
> gauss_base_func <- function(x) exp( -(x - u_i)^2/(2*s*s) )
> plot(gauss_base_func)

f:id:n_shuyo:20190427122025p:plain

「空間的な尺度」とやらを表すパラメータ s について、PRML 本では特に指定がない。
ちなみにコンソールで対話的に定義した変数はグローバル。s を変えて plot() し直すとグラフに反映されるので、塩梅のいい s を探す。
0.1 でちょうど本に載ってるグラフと同じような形になったので、これでいいや。


「空間的な尺度」を表すパラメータ μ_i を [0, 1] で９等分にとって、ガウス基底関数を９つ用意する必要がある。
ここで「スカラー値をとる関数にベクトルを与えると、各要素に関数を適応したベクトルを返す」というＲの特徴を思い出すと、こんな風に書くこともできる。

> phi <- function(x)exp(-(x-seq(0,1,length=9))^2/(2*s*s))
> phi(0)
[1] 1.000000e+00 4.578334e-01 4.393693e-02 8.838263e-04 3.726653e-06
[6] 3.293714e-09 6.101937e-13 2.369542e-17 1.928750e-22


u_i のところを seq(0,1,length=9)、[0, 1] で９等分した点から成る「ベクトル」に変えると、スカラー値を与えるとベクトルが返る関数が書けてしまう。
Ｒバンザイ。


ではこれに n次のベクトルを与えると、 n×9 の行列が！　というわけにはいかない。

> phi(c(0, 0.5))
[1] 1.000000e+00 8.838263e-04 4.393693e-02 4.578334e-01 3.726653e-06
[6] 4.578334e-01 6.101937e-13 8.838263e-04 1.928750e-22
Warning message:
In x - seq(0, 1, length = 9) :
   長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません


またしばらく考え＆試行錯誤。

> sapply(c(1, 0.5), phi)
              [,1]         [,2]
 [1,] 1.928750e-22 3.726653e-06
 [2,] 2.369542e-17 8.838263e-04
 [3,] 6.101937e-13 4.393693e-02
 [4,] 3.293714e-09 4.578334e-01
 [5,] 3.726653e-06 1.000000e+00
 [6,] 8.838263e-04 4.578334e-01
 [7,] 4.393693e-02 4.393693e-02
 [8,] 4.578334e-01 8.838263e-04
 [9,] 1.000000e+00 3.726653e-06


うおおおおおお、できた。
一目瞭然とまではいかないけど、でも数学の数式と同じ行数で design matrix Φ を求められるなんて。
Ｒすごいな。


phi() にフィードバック。
今後はベクトルを突っ込んでも大丈夫。

> phi <- function(x) sapply(x, function(x){ exp(-(x-seq(0,1,length=9))^2/(2*s*s)) })

ややこしすぎて読めないって？　見慣れたらそうでもないｗ
まあでも sapply 無い方がよっぽど見通しいいよなあ。
match.fun とか駆使したらもっとうまく書けたりするんだろうか。



ここまでくれば線形回帰はできたも同然だけど、データがないから作らなきゃ。
データフレームを使うのが清く正しいっぽい。


PRML 本には \sin\left(2\pi x\right) にガウスノイズを加え、と例によって書いてあるが、ノイズの分散は指定されていない。
3.3.1 と同じ 0.2 にしておくか。他のハイパーパラメータα、βも同じくそれぞれ 2, 25 としておく。

> xlist <- seq(0, 1, length=25)
> tlist <- sin(2*pi*xlist)+rnorm(length(xlist), sd=0.2)
> D0 <- data.frame(x=xlist, t=tlist)
> D0
            x           t
1  0.00000000 -0.06303114
2  0.04166667 -0.18576933
3  0.08333333  0.68159289
        :           :
23 0.91666667 -0.67836209
24 0.95833333 -0.44500192
25 1.00000000  0.17576182
> plot(D0)

f:id:n_shuyo:20190427122052p:plain


ここでちょっとした罠が。
"sin(2*pi*xlist)+rnorm(1, sd=0.2)" と書いてしまうと、rnorm(1, sd=0.2) で得られた「１個の」ノイズ値が２５個コピーされて加算されてしまう。
"rnorm(length(xlist), sd=0.2)" と書かないといけない。


標本をランダムに選ばせるには sample() が使える。

> (D <- D0[sample(nrow(D0), 4),])  # choice 4 training data at random
            x           t
8  0.29166667  1.01907177
14 0.54166667 -0.04237673
2  0.04166667 -0.18576933
21 0.83333333 -0.76052558

データもできたので、いよいよ予測分布を計算してみる。
ここまで揃っていれば後は簡単、なはず。


計画行列(design matrix) Φ を計算。

> phi <- function(x)sapply(x,function(x)exp(-(x-seq(0,1,length=9))^2/(2*s*s)))
> PHI <- phi(D$x)
> PHI
             [,1]         [,2]         [,3]        [,4]        [,5]         [,6]         [,7]         [,8]         [,9]
[1,] 1.421479e-02 2.493522e-01 9.168554e-01 0.706648278 0.114161760 3.865920e-03 2.744100e-05 4.082836e-08 1.273324e-11
[2,] 4.254413e-07 1.698567e-04 1.421479e-02 0.249352209 0.916855356 7.066483e-01 1.141618e-01 3.865920e-03 2.744100e-05
[3,] 9.168554e-01 7.066483e-01 1.141618e-01 0.003865920 0.000027441 4.082836e-08 1.273324e-11 8.323970e-16 1.140610e-20
[4,] 8.323970e-16 1.273324e-11 4.082836e-08 0.000027441 0.003865920 1.141618e-01 7.066483e-01 9.168554e-01 2.493522e-01


事後分布の期待値 m_N、共分散 S_N を計算。

> S_N_inv <- alpha * diag(9) + beta * t(PHI) %*% PHI
> S_N <- solve(S_N_inv)
> m_N <- beta * S_N %*% t(PHI) %*% D$t


回帰関数を求めて、グラフ。
回帰関数が簡潔に求められているのが、なかなか気持ちいい。

> y <- function(x) (t(phi(x)) %*% m_N)             # regression function
> plot(y, xlim=c(0,1), ylim=c(-1.2,1.2))
> par(new=T)
> plot(D, xlim=c(0,1), ylim=c(-1.2,1.2), ylab="")  # distribution of data


予測分布を求めて、グラフ。
mapply は無くてもうまくいきそうな気もするのだけど、今のところ妙案を思いつかず。

var_N <- function(x){ 1/beta + (t(phi(x)) %*% S_N %*% phi(x))[1] }

# predictive distribution
p <- function(x, t){
  mapply(function(x, t)dnorm(t, m=(t(m_N) %*% phi(x))[1], s=var_N(x), log=T), x, t)
}

x <- seq(0,1,length=30)
t <- seq(-1.5, 1.5, length=60)
z <- outer(x, t, p)
persp(x, t, z, theta=90, phi=60, shade=0.4)


ベクトルの内積が (t(v1) %*% v2)[1] となっているのがダサい。他にいい書き方ないかなあ。
【追記】 ああ sum(v1*v2) でもいいのか。これもいまいち っちゃあ いまいちだけど、タイプ数はだいぶ少なくなった。



全部つなげて完成。
対話的にいろいろな条件で試しやすいよう、２つの関数にまとめる。


グラフ描画とかあれこれで長くなっているが、本質的な部分は 10 行程度。
でもまだもうちょっと短くできそうだな……

M <- 9        # number of basis functions
alpha <- 2    # hyper parameter
beta <- 25    # hyper parameter(precision)
Lattice <- 30 # number of graph's lattice
s <- 0.1

# training data
xlist <- seq(0, 1, length=25)
tlist <- sin(2*pi*xlist)+rnorm(length(xlist), sd=0.2)
D0 <- data.frame(x=xlist, t=tlist)

predictive <- function(D) {
    # design matrix
    phi <- function(x) sapply(x,function(x)exp(-(x-seq(0,1,length=9))^2/(2*s*s)))
    PHI <- t(phi(D$x))

    # covariance matrix & means
    S_N_inv <- alpha * diag(M) + beta * t(PHI) %*% PHI
    S_N <- solve(S_N_inv)
    m_N <- beta * S_N %*% t(PHI) %*% D$t

    # regression function
    y <- function(x) (t(phi(x)) %*% m_N)
    plot(y, xlim=c(0,1), ylim=c(-1.2,1.2))
    par(new=T)
    plot(D, xlim=c(0,1), ylim=c(-1.2,1.2), ylab="")

    # predictive distribution
    var_N <- function(x){ 1/beta + (t(phi(x)) %*% S_N %*% phi(x))[1] }
    function(x, t){ 
        mapply(function(x, t)dnorm(t, m=(t(m_N) %*% phi(x))[1], s=var_N(x), log=T), x, t)
    }
}
draw_dist <- function(p){
    x <- seq(0,1,length=Lattice)
    t <- seq(-1.5, 1.5, length=Lattice*2)
    z <- outer(x, t, p)
    persp(x, t, z, theta=0, phi=60, shade=0.4)
}

p <- predictive(D0[sample(length(D0$x))[1:1],])
draw_dist(p);
p <- predictive(D0[sample(length(D0$x))[1:2],])
draw_dist(p);
p <- predictive(D0[sample(length(D0$x))[1:4],])
draw_dist(p);
p <- predictive(D0)
draw_dist(p);

f:id:n_shuyo:20190427122141p:plainf:id:n_shuyo:20190427122149p:plainf:id:n_shuyo:20190427122158p:plainf:id:n_shuyo:20190427122207p:plain

訓練データがあるところは分布が削れてるんだけど、予測分布の分散には実はカーネル関数が現れてるから……とかいろいろ ふむふむポイントが多くて、楽しい。


Ｒは……まだ全然使いこなせている実感はないけど、ちょっと試してみたいこと(PRML読んでりゃいっぱいあるでしょー)を気軽に試せる手段ができて、いいね。