Nothing
index.map <- function(L, l, k){
return((2*L-k)*(k-1)/2+l)
}
relevant.funs <- function(intercept=TRUE){
## high dimensional linear regression ##
train.fun <- function(X, y, lambda=NULL){
if(is.null(lambda)) lambda = 'CV.min'
p = ncol(X)
htheta <- if (lambda == "CV.min") {
outLas <- cv.glmnet(X, y, family = "gaussian", alpha = 1,
intercept = intercept, standardize = T)
as.vector(coef(outLas, s = outLas$lambda.min))
} else if (lambda == "CV") {
outLas <- cv.glmnet(X, y, family = "gaussian", alpha = 1,
intercept = intercept, standardize = T)
as.vector(coef(outLas, s = outLas$lambda.1se))
} else {
outLas <- glmnet(X, y, family = "gaussian", alpha = 1,
intercept = intercept, standardize = T)
as.vector(coef(outLas, s = lambda))
}
if(intercept==FALSE) htheta = htheta[2:(p+1)]
return(list(lasso.est = htheta))
}
## predict ##
pred.fun <- function(X, htheta){
X = as.matrix(X)
if((length(htheta) - ncol(X)) == 1) X = cbind(1,X)
pred = as.vector(X%*%htheta)
return(pred)
}
## variance of residual ##
dev.fun <- function(pred, y, sparsity=0){
pred = as.vector(pred); y = as.vector(y)
n = length(y)
sigmasq.hat = sum((y - pred)^2) / max(0.7*n, n-sparsity)
return(sigmasq.hat)
}
return(list(train.fun = train.fun,
pred.fun = pred.fun,
dev.fun = dev.fun))
}
opt.weight<-function(Gamma,delta,report.reward=TRUE){
## Purpose: Compute Ridge-type weight vector
## Returns: weight: the minimizer \eqn{\gamma}
## reward: the value of penalized reward
## ----------------------------------------------------
## Arguments: Gamma: regression covariance matrix, of dimension \eqn{L} x \eqn{L}
## delta the ridge penalty level, non-positive.
## report.reward the reward is computed or not (Default = `TRUE`)
## ----------------------------------------------------
L<-dim(Gamma)[2]
opt.weight<-rep(NA, L)
opt.reward<-NA
# Problem definition
v<-Variable(L)
Diag.matrix<-diag(eigen(Gamma)$values)
for(ind in 1:L){
Diag.matrix[ind,ind]<-max(Diag.matrix[ind,ind],0.001)
}
Gamma.positive<-eigen(Gamma)$vectors%*%Diag.matrix%*%t(eigen(Gamma)$vectors)
objective <- Minimize(quad_form(v,Gamma.positive+diag(delta,L)))
constraints <- list(v >= 0, sum(v)== 1)
prob.weight<- Problem(objective, constraints)
if(is_dcp(prob.weight)){
result<- solve(prob.weight)
opt.status<-result$status
opt.sol<-result$getValue(v)
for(l in 1:L){
opt.weight[l]<-opt.sol[l]*(abs(opt.sol[l])>10^{-8})
}
}
if(report.reward){
v<-Variable(L)
objective<-Minimize(2*t(v)%*%Gamma.positive%*%opt.weight-t(opt.weight)%*%Gamma.positive%*%opt.weight)
constraints<-list(v >= 0, sum(v)== 1)
delta.optim<-Problem(objective, constraints)
result<- solve(delta.optim)
opt.reward<-result$value
returnList <- list("weight" = opt.weight,
"reward" = opt.reward)
}else{
returnList <- list("weight" = opt.weight)
}
return(returnList)
}
gensamples <- function(gen.mu, gen.Cov, gen.size=500, threshold=0, alpha.thres=0.01){
gen.dim = length(gen.mu)
if(threshold==0){
thres = qnorm(1-alpha.thres/(gen.dim*2))
gen.samples = matrix(0, nrow=gen.size, ncol=gen.dim)
n.picked = 0 # records how many samples are picked
while(n.picked < gen.size){
S = mvrnorm(1, mu=rep(0, gen.dim), Sigma=gen.Cov)
if(max(abs(S / sqrt(diag(gen.Cov)))) <= thres){
n.picked = n.picked + 1
gen.samples[n.picked, ] = gen.mu + S
}
}
}
if(threshold==1){
# 1 stands for chi square threshold
UDV_list = svd(gen.Cov)
U = UDV_list$u
D = UDV_list$d
V = UDV_list$v
D.sqrt = sqrt(D)
gen.Cov.sqrt = U %*% diag(D.sqrt) %*% t(V)
thres = 1.0 * qchisq((1-alpha.thres), df=gen.dim)
gen.samples = matrix(0, nrow=gen.size, ncol=gen.dim)
n.picked = 0
while(n.picked < gen.size){
Z = mvrnorm(1, mu=rep(0, gen.dim), Sigma=diag(gen.dim))
Z.normsq = sum(Z^2)
if(Z.normsq <= thres){
n.picked = n.picked + 1
gen.samples[n.picked, ] = gen.mu + gen.Cov.sqrt %*% Z
}
}
}
if(threshold==2){
gen.samples = matrix(mvrnorm(gen.size, gen.mu, Sigma=gen.Cov),
nrow=gen.size, ncol=gen.dim)
}
return(gen.samples)
}
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.