Nothing
picasso.sqrtlasso <- function(X,
Y,
lambda = NULL,
nlambda = NULL,
lambda.min.ratio = NULL,
method="l1",
gamma = 3,
standardize = TRUE,
intercept = TRUE,
prec = 1e-4,
max.ite = 1e4,
verbose = FALSE)
{
n = nrow(X)
d = ncol(X)
Yb = Y
begt = Sys.time()
if (verbose)
cat("Sparse logistic regression. \n")
if (n == 0 || d == 0) {
cat("No data input.\n")
return(NULL)
}
if (method!="l1" && method!="mcp" && method!="scad"){
cat(" Wrong \"method\" input. \n \"method\"
should be one of \"l1\", \"mcp\", \"scad\".\n",
method,"is not supported in this version. \n")
return(NULL)
}
if (standardize){
xx = rep(0.0, n*d)
xm = rep(0.0, d)
xinvc.vec = rep(0.0, d)
str = .C("standardize_design", as.double(X), as.double(xx), as.double(xm), as.double(xinvc.vec),
as.integer(n), as.integer(d), PACKAGE="picasso")
xx = matrix(unlist(str[2]), nrow=n, ncol=d, byrow=FALSE)
xm = matrix(unlist(str[3]), nrow=1)
xinvc.vec = unlist(str[4])
}else{
xinvc.vec = rep(1,d)
xx = X
}
yy = Yb
if (!is.null(lambda)) nlambda = length(lambda)
if (is.null(lambda)){
if (is.null(nlambda))
nlambda = 100
L0 = sqrt(sum(yy*yy)/n)
lambda.max = max(abs(crossprod(xx,yy/n)))/L0
if (is.null(lambda.min.ratio)){
lambda.min = 0.05*lambda.max
} else {
lambda.min = min(lambda.min.ratio*lambda.max, lambda.max)
}
if (lambda.min >= lambda.max)
cat("lambda.min is too small. \n")
lambda = exp(seq(log(lambda.max), log(lambda.min), length = nlambda))
}
if (method == "l1" || method == "mcp" || method == "scad") {
if(method == "l1") {
method.flag = 1
}
if (method == "scad") {
method.flag = 3
if (gamma <= 2) {
cat("gamma > 2 is required for SCAD. Set to default value 3. \n")
gamma = 3
}
}
if (method == "mcp") {
method.flag = 2
if (gamma <= 1) {
cat("gamma > 1 is required for MCP. Set to default value 3. \n")
gamma = 3
}
}
out = sqrtlasso_solver(yy, xx, lambda, nlambda, gamma,
n, d, max.ite, prec, intercept, verbose,
method.flag)
}
df = rep(0, nlambda)
for (i in 1:nlambda)
df[i] = sum(out$beta[[i]]!=0)
est = list()
intcpt = matrix(0, nrow=1, ncol=nlambda)
beta1 = matrix(0, nrow=d, ncol=nlambda)
if (standardize){
for (k in 1:nlambda){
tmp.beta = out$beta[[k]]
beta1[,k] = xinvc.vec*tmp.beta
intcpt[k] = -as.numeric(xm[1,]%*%beta1[,k])+out$intcpt[k]
}
} else {
for (k in 1:nlambda){
beta1[,k] = out$beta[[k]]
intcpt[k] = out$intcpt[k]
}
}
runt = Sys.time()-begt
est$runt = out$runt
est$beta = Matrix(beta1)
res = X%*%beta1+matrix(rep(intcpt,n),nrow=n,byrow=TRUE)
est$p = exp(res)/(1+exp(res))
est$intercept = intcpt
est$lambda = lambda
est$nlambda = nlambda
est$df = df
est$method = method
est$ite =out$ite
est$verbose = verbose
est$runtime = runt
class(est) = ""
return(est)
}
print.sqrtlasso <- function(x, ...)
{
cat("\n SQRT Lasso options summary: \n")
cat(x$nlambda, " lambdas used:\n")
print(signif(x$lambda,digits=3))
cat("Degree of freedom:",min(x$df),"----->",max(x$df),"\n")
if(units.difftime(x$runtime)=="secs") unit="secs"
if(units.difftime(x$runtime)=="mins") unit="mins"
if(units.difftime(x$runtime)=="hours") unit="hours"
cat("Runtime:",x$runtime," ",unit,"\n")
}
plot.sqrtlasso <- function(x, ...)
{
matplot(x$lambda, t(x$beta), type="l", main="Regularization Path",
xlab="Regularization Parameter", ylab="Coefficient")
}
coef.sqrtlasso <- function(object, lambda.idx = c(1:3), beta.idx = c(1:3), ...)
{
lambda.n = length(lambda.idx)
beta.n = length(beta.idx)
cat("\n Values of estimated coefficients: \n")
cat(" index ")
for(i in 1:lambda.n){
cat("",formatC(lambda.idx[i],digits=5,width=10),"")
}
cat("\n")
cat(" lambda ")
for(i in 1:lambda.n){
cat("",formatC(object$lambda[lambda.idx[i]],digits=4,width=10),"")
}
cat("\n")
cat(" intercept ")
for(i in 1:lambda.n){
cat("",formatC(object$intercept[i],digits=4,width=10),"")
}
cat("\n")
for(i in 1:beta.n){
cat(" beta",formatC(beta.idx[i],digits=5,width=-5))
for(j in 1:lambda.n){
cat("",formatC(object$beta[beta.idx[i],lambda.idx[j]],digits=4,width=10),"")
}
cat("\n")
}
}
predict.sqrtlasso <- function(object, newdata, lambda.idx = c(1:3), Y.pred.idx = c(1:5), ...)
{
pred.n = nrow(newdata)
lambda.n = length(lambda.idx)
Y.pred.n = length(Y.pred.idx)
intcpt = matrix(rep(object$intercept[lambda.idx],pred.n),nrow=pred.n,
ncol=lambda.n,byrow=T)
Y.pred = newdata%*%object$beta[,lambda.idx] + intcpt
cat("\n Values of predicted responses: \n")
cat(" index ")
for (i in 1:lambda.n){
cat("",formatC(lambda.idx[i], digits=5, width=10),"")
}
cat("\n")
cat(" lambda ")
for (i in 1:lambda.n){
cat("",formatC(object$lambda[lambda.idx[i]], digits=4, width=10),"")
}
cat("\n")
for (i in 1:Y.pred.n){
cat(" Y",formatC(Y.pred.idx[i], digits=5, width=-5))
for (j in 1:lambda.n){
cat("",formatC(Y.pred[Y.pred.idx[i],j], digits=4, width=10),"")
}
cat("\n")
}
return(Y.pred)
}
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.