intro_to_ml/23_exercices.R

85 lines
1.3 KiB
R

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
}