# set-up the design matrix
p <- 21; n <- (p-1)*5+1; x <- matrix(0,n,p)
for (i in 1:(p-1)) {
 x[(5*i-4):(5*i),i] <- seq(from=1,to=.2,by=-.2)
 x[(5*i-4):(5*i),i+1] <- seq(from=0,to=.8,by=.2) }
x[n,p] <- 1

# set-up the penalty matrix (hope I got the entries right!)
mu <- rep(0,p)
b <- matrix(0,p,p)
for (i in 3:(p-2)) {
  b[i,(i-2):(i+2)] <- c(.25,-1,1.5,-1,.25) }
b[1,1:3] <- c(.25, -.5, .25); 
b[2,1:4] <- c(-.5,1.25,-1,.25)
b[p,(p-2):p] <- c(.25,-.5,.25)
b[p-1,(p-3):p] <- c(.25,-1,1.25,-.5)

# set-up "true" parameter values and generate data
theta <- rep(0,21)
theta[1:5] <- 2*(1:5)
theta[6:10] <- 10+.8*(1:5)
theta[11:15] <- 14+.2*(1:5)
theta[16:21] <- 15
set.seed(133)
y <- rnorm(n, x%*%theta, 5)


ps.options(horizontal=T)
postscript("penalty.ps")
par(mfrow=c(3,3))

# plot data
plot((0:(n-1))*((p-1)/(n-1)), y, xlab="t",ylab="y") 

### try different strengths of penalty
for (lambda in c(0, 5, 10, 30, 50, 200, 2000, 5000)) {

### remember method 2 from last time to compute fit
u <- chol( t(x) %*% x + lambda*b)  
tmp <-  backsolve(u, t(x) %*% y + lambda*b%*%mu, transpose=T)
ans <-  backsolve(u, tmp, transpose=F) 

### how much wiggle in fit?
print(c(lambda,
        sum( (ans[2:(p-1)]-.5*ans[1:(p-2)]-.5*ans[3:p])^2 ) ), digits=3)

### plot estimated f()
plot(0:(p-1), ans, type="l", ylim=range(y),xlab="t",ylab="y",col=4)
title(paste("lambda =",as.character(lambda)))
### superimpose true f()
points(0:(p-1), theta, type="l", lty=2,col=2)

}

graphics.off()