R 정규화 최소 2승법

PRML3.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)