Nothing
### functions that are called by the main functions, but invisible to users unless using namespace ':::'
## standardize the data ------------------------------------------------------------
# standardize x to be mean 0 and norm 1
# standardize y to be mean 0
std <- function(x, y, intercept){
n = dim(x)[1]
p = dim(x)[2]
if(intercept){
mean_x = Matrix::colMeans(x)
mean_y = mean(y)
x = scale(x, center = mean_x, scale = FALSE)
y = scale(y, center = mean_y, scale = FALSE)
}else{
mean_x = rep(0, p)
mean_y = 0
}
sd_demeanedx = sqrt(Matrix::colSums(x^2))
x = scale(x, center = FALSE, scale = sd_demeanedx)
return(list(x_std = x, y_std = y, mean_x=mean_x, sd_demeanedx=sd_demeanedx, mean_y=mean_y))
}
## calculate various information criteria: AIC, BIC, AICc, BICc, GCV, Cp -----------
calc.ic.all <- function(coef, x, y, df, sigma=NULL){
# unify dimensions
y = matrix(y, ncol=1)
df = matrix(df, nrow=1)
# sanity check
if(ncol(coef) != ncol(df)){
stop('the number of coef vectors does not match the number of df')
}
if(is.null(sigma)){
stop("need to specify sigma for Mallow's Cp")
}
n = nrow(y)
nfit = ncol(coef)
fit = x %*% coef
rss = Matrix::colSums(sweep(fit, 1, y, '-')^2)
ic = list()
ic$aic = log(rss/n) + 2*df/n
ic$bic = log(rss/n) + log(n)*df/n
ic$gcv = rss / (n-df)^2
ic$cp = rss + 2*sigma^2*df
# for AICc and BICc df larger than n-2 will cause trouble, round it
df[which(df>=n-2)] = n-3
ic$aicc = log(rss/n) + 2*(df+1)/(n-df-2)
ic$bicc = log(rss/n) + log(n)*(df+1)/(n-df-2)
# ic$bicc = log(rss/n) + df*(log(n)+2*log(p))/n
return(ic)
}
## Heuristic df for BOSS ------------------------------------------------------------
# @param Q An orthogonal matrix, with \code{nrow(Q)=length(y)=n} and
# \code{ncol(Q)=p}. For BOSS, Q is obtained by QR decomposition upon
# an ordered design matrix.
# @param y A vector of response variable, with \code{length(y)=n}.
# @param sigma,mu The standard deviation and mean vector of the true model.
# In practice, if not specified, they are calculated via full multiple
# regression of y upon Q.
# @param x is required for the case where n<=p and sigma is estimated via lasso estimate
calc.hdf <- function(Q, y, sigma=NULL, mu=NULL, x=NULL){
n = dim(Q)[1]
p = dim(Q)[2]
# if(p>=n){
# stop('hdf is undefined when p>=n')
# }
if(n <= p & is.null(x)){
stop("This is the n<p scenario. Need the design matrix X for the lasso estimated sigmahat.")
}
# mu is estimated via multiple LS regression
beta_hat = t(Q) %*% y # the multiple regression coef
if(is.null(mu)){
xtmu = beta_hat
}else{
xtmu = t(Q)%*%mu
}
xtmu_matrix = matrix(rep(xtmu,each=p-1), ncol=p-1, byrow=TRUE)
if(is.null(sigma)){
# For n>p, the estimated sigma is based on multiple LS regression
# For n<=p, the estimated sigma is based on lasso with 10-fold CV
if(n > p){
resid = y - Q %*% beta_hat
sigma = sqrt(sum(resid^2)/(n-p))
}else{
sigma = est.sigma.lasso(x, y)
}
}
tryCatch({
# calculate the inverse function of E(k(lambda))=k, where k=1,...p-1
inverse = function(f, lower, upper) {
function(y) stats::uniroot(function(x){f(x) - y}, lower=lower, upper=upper)[1]
}
exp_size <- function(x){
c = stats::pnorm((x-xtmu) / sigma)
d = stats::pnorm((-x-xtmu) / sigma)
return( sum(1 - c + d) )
}
inverse_exp_size = inverse(exp_size, 0, 100*max(abs(xtmu)))
sqrt_2lambda = unlist(lapply(1:(p-1), inverse_exp_size))
sqrt_2lambda_matrix = matrix(rep(sqrt_2lambda,each=p), nrow=p, byrow=F)
# plug the sequence of lambda into the expression of df(lambda)
a = stats::dnorm((sqrt_2lambda_matrix-xtmu_matrix) / sigma)
b = stats::dnorm((-sqrt_2lambda_matrix-xtmu_matrix) / sigma)
size = 1:(p-1)
sdf = (sqrt_2lambda/sigma) * Matrix::colSums(a + b)
df = size + sdf
names(df) = NULL
return(list(hdf=c(0, df, p), sigma=sigma))
}, error=function(e){
warning('returns the df for a null model')
return(list(hdf=c(0, 1:p + 2*p*stats::qnorm(1-1:p/(2*p))*stats::dnorm(stats::qnorm(1-1:p/(2*p)))), sigma=sigma))
})
}
## Estimate sigma based on lasso-cv
est.sigma.lasso <- function(x, y){
n = dim(x)[1]
lasso_cv = glmnet::cv.glmnet(x, y, intercept=FALSE)
betahat = glmnet::coef.glmnet(lasso_cv, s='lambda.min')[-1,]
df = sum(betahat != 0)
y_hat = x %*% betahat
sigma_hat = sqrt( sum((y-y_hat)^2) / (n-df-1) )
return(sigma_hat)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.