rm(list=ls())
set.seed(1123)

source('15_loocv.R')

multdiag <-
function(X,d)
{
  R <- matrix(NA, nrow=dim(X)[1], ncol=dim(X)[2])
  for (i in 1:dim(X)[2]) { R[,i]=X[,i]*d[i] }
  return(R)
}

n <- 700
p <- 55

sd <- 6 # standard deviation for zero-mean gaussian noise
X <- matrix(runif(n*p),nrow=n,ncol=p)
X <- scale(X)
beta <- runif(p, min=-10, max=10)
y <- X%*%beta + rnorm(n, mean=0, sd=sd)

lambdas <- 10^seq(-1,3,by=0.2)

var <-
function(lambda)
{
  d <- (Xs$d^2)/(Xs$d^2 + lambda)^2
  var <- multdiag(Xs$v,d)
  var <- sd^2 * tcrossprod(var,Xs$v)
}

bias <-
function(lambda)
{
  d <- lambda/(Xs$d^2+lambda)
  bias <- multdiag(Xs$v,d)
  bias <- bias %*% crossprod(Xs$v,beta)
}

epeVar <-
function(lambda)
{
  var <- var(lambda)
  return( mean(rowSums(X*(X%*%var))) )
}

epeBias <-
function(lambda)
{
  bias <- bias(lambda)
  return( mean((X%*%bias)^2) )
}

epe <-
function(lambda)
{
  return( epeVar(lambda) + epeBias(lambda) + sd^2 )
}

#rm <- ridge(X, y, lambdas)
#Xs <- svd(X)
#epes <- sapply(lambdas, epe)

X.init <- X
beta.init <- beta

ps <- seq(1,p)
lps <- length(ps)
epes <- numeric(lps)
biass <- numeric(lps)
vars <- numeric(lps)
maes <- numeric(lps)

for (k in ps)
{
  X <- X.init[,1:k]
  beta <- beta.init[1:k]
  rm <- ridge(X, y, lambdas)
  Xs <- svd(X)
  epes[k] <- epe(rm$lambda)
  biass[k] <- epeBias(rm$lambda)
  vars[k] <- epeVar(rm$lambda)
  maes[k] <- rm$mae
}