85 lines
1.3 KiB
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
|
|
} |