Rで正則化最小二乗法

PRML 3.1.4に記載されているように、重みベクトルの二乗和を正則化項とした誤差関数を最小化するような、線形基底関数モデル（多項式基底）の係数を求め、正則化パラメータの値によってモデルの過学習が抑制される様子を図にします。

set.seed(0)
x<-seq(0.05,1,0.1)
t<-sin(2*pi*x)+rnorm(10,0,0.5)
t<-matrix(t,nrow=length(t))
M<-10
base<-function(m,x) {
        r<-x^m                          #多項式基底の場合
#        r<-exp(-(x-(m-1)/(M-2))^2/(2*0.2^2))    #ガウス基底の場合
        r[m == 0]<-1                    # w0をバイアス項とするためφ0=1
        r
}
makePhi<-function(M) {
        A <- matrix(nrow=length(x),ncol=M)
        for (i in 1:length(x)) {
                for (j in 0:(M-1)) {
                        A[i,1+j] <- base(j,x[i])
                }
        }
        A
}
xrange<-c(-0.5,1.5)
yrange<-c(-2,2)
plot(x,t,xlim=xrange,ylim=yrange)
phi<-makePhi(M)
lambdabase<-10:0
for (ll in lambdabase) {
        lambda<-10^(-ll*2)
        w<-solve(lambda * diag(M) + t(phi) %*% phi) %*% t(phi) %*% t
        print(w)
        estimate<-function(x){
                i <- 0:(M-1)
                (sapply(x,function(xn){ (sum(w * base(i, xn))) } )) # (w[1]+w[2]*base(1,x)+w[3]*base(2,x)...)
        }
        par(new=T)
        curve(estimate,type="l",xlim=xrange,ylim=yrange,col=ll+1)
}
legend("topleft",legend=10^(-lambdabase*2),col=lambdabase+1,lty=1,cex=0.6)