線形回帰モデル

PRML読書会の予習を兼ねて3.1のガウス基底関数を用いて線形回帰を行うプログラムをRで書いた。

疑問としては(3.4)式の

の部分で\mu_jとsはどうやって決定するのかがよくわからなかった。とりあえずs=1,\mu_j = j / Mとした。

追記：本当はbase(x[i],j)ではなくてbase[ [j] ](x[i])と書きたいんだけどRで

for(j in 1:10){
  base[[j]] <- function(x) return x^i
}

とか書くと全部10乗を返す関数になってしまうのでbase(x[i],j)とか書いてるけどなんかいい方法はないものか。

#generate data set
gen_data <- function(N = 25){
  gauss_noise <- rnorm(N)
  x <- runif(N)
  y <- numeric(N)
  sincurve<-function(x) sin(2 * pi *x)
  for(i in 1:N){
    y[i] <- sincurve(x[i]) + 0.3 * gauss_noise[i]
  }
  return (list(x,y))
}

#x: input value
#y: target value
#M: order
#base: base function
#lambda: regularization parameter
lreg <- function(x , y , M , base , lambda = 0){
  N <- length(x)
  I <- matrix(0 , M , M)
  dmat <- matrix(0 , N , M)
  for(i in 1:N){
    for(j in 1:M){
      dmat[i,j] = base(x[i], j)
    }
  }
  inv <- solve(t(dmat) %*% dmat + lambda * diag(M))
  w <- inv %*% (t(dmat) %*% y)
  reg.function <- function(x){
    res = 0.0
    for(i in 1:M){
      res =res + w[i] * base(x , i)
    }
    return (res)
  }
  return (list(w , reg.function))
}

poly_base <- function(x, i){
  return (x^(i-1))
}

l <- gen_data(25)
x <- l[[1]]
y <- l[[2]]
M = 25
gauss_base <- function(x,i){
  if(i == 1){
    return (0)
  }
  return (exp(-((x-1.0 * i / M)^2)))
}
ws <- lreg(x,y,25,gauss_base , 0.0001)
plot(x,y)
plot(ws[[2]] , add=T)