### SVM based on composite operator irsvm <- function(x, ...) UseMethod("irsvm") irsvm.default <- function(x, ...) { if (extends(class(x), "Matrix")) return(irsvm.matrix(x = x, ...)) if (inherits(x, "data.frame")) return(irsvm.matrix(x = x, ...)) if (inherits(x, "numeric") && is.null(dim(x))) return(irsvm.matrix(x = as.matrix(x), ...)) stop("no method for objects of class ", sQuote(class(x)), " implemented") } irsvm.formula <- function(formula, data, weights, contrasts=NULL, ...){ ## extract x, y, etc from the model formula and frame if(missing(data)) data <- environment(formula) mf <- match.call(expand.dots = FALSE) m <- match(c("formula", "data", "weights"), names(mf), 0L) mf <- mf[c(1L, m)] mf$drop.unused.levels <- TRUE mf[[1L]] <- as.name("model.frame") mf <- eval(mf, parent.frame()) mt <- attr(mf, "terms") # allow model.frame to have updated it Y <- model.response(mf, "any") # e.g. factors are allowed ## avoid problems with 1D arrays, but keep names if(length(dim(Y)) == 1L) { nm <- rownames(Y) dim(Y) <- NULL if(!is.null(nm)) names(Y) <- nm } ## null model support X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts) else matrix(,NROW(Y), 0L) ## avoid any problems with 1D or nx1 arrays by as.vector. weights <- as.vector(model.weights(mf)) if(!length(weights)) weights <- rep(1, nrow(mf)) else if(any(weights < 0)) stop("negative weights not allowed") if(!is.null(weights) && !is.numeric(weights)) stop("'weights' must be a numeric vector") if(length(weights) != length(Y)) stop("'weights' must be the same length as response variable") RET <- irsvm_fit(X, Y, weights,...) RET$call <- match.call() RET <- c(RET, list(formula=formula, terms = mt, data=data, contrasts = attr(X, "contrasts"), xlevels = .getXlevels(mt, mf))) class(RET) <- "wsvm" RET } irsvm.matrix <- function(x, y, weights, ...){ RET <- irsvm_fit(x, y, weights,...) RET$call <- match.call() return(RET) } irsvm_fit <- function(x,y, weights, cfun="ccave", s=NULL, delta=0.0001, type=NULL, kernel = "radial", #degree = 3, #gamma = if (is.vector(x)) 1 else 1 / ncol(x), #coef0 = 0, cost = 1, #class.weights = NULL, #cachesize = 100, #tolerance = 0.001, epsilon = 0.1, #shrinking = TRUE, #probability = FALSE, #fitted = TRUE, #subset, #na.action, iter=10, reltol=1e-5, trace=FALSE, ...){ call <- match.call() ## determine model type if (is.null(type)) type <- if (is.factor(y)) "C-classification" else "eps-regression" cfunold <- cfun cfun <- cfun2num(cfun) dfun <- switch(type, "C-classification"=6, "nu-classification"=6, "eps-regression"=2, "nu-regression"=2) if(dfun==6){ if (is.factor(y)) { y <- as.integer(y) } else { if(any(as.integer(y) != y)) stop("y variable has to be of factor or integer type for classification\n") y <- as.factor(y) y <- as.integer(y) } if(length(unique(y))!=2) stop("y must be a binary variable for classification\n") else if(!all(names(table(y)) %in% c(1, -1))){ y[y==1] <- -1 y[y==2] <- 1 } } if(is.null(s)) s <- assign_s(cfun, y) else check_s(cfun, dfun=6, s) if(cfun==6) if(s > 1) delta <- (s-1)/2 else if(s==1) delta <- 0 else{ if(is.null(delta)) stop("delta must be provided") if(delta <= 0) stop("delta must be positive") } ## avoid any problems with 1D or nx1 arrays by as.vector. if(missing(weights)) weights=rep(1,length(y)) weights <- as.vector(weights) if(!is.null(weights) && !is.numeric(weights)) stop("'weights' must be a numeric vector") ## check weights if( !is.null(weights) && any(weights < 0) ){ stop("negative weights not allowed") } if(length(weights) != length(y)) stop("'weights' must be the same length as response variable") n <- length(y) if(all(x[,1]==1)) xtmp <- x[,-1] else xtmp <- x d1 <- 10 k <- 1 if(trace) { cat("\niterations ...\n") } los <- rep(NA, iter) weights_update <- weights subset <- 1:n while(d1 > reltol && k <= iter){ RET <- wsvm(xtmp, y, weights*weights_update, type=type, kernel=kernel, cost=cost, epsilon=epsilon, cross=0, subset=subset, ...) pred <- predict(RET, xtmp, decision.values=TRUE) if(dfun==6){ fitted.values <- attr(pred, "decision.values") if(y[subset[1]]==-1) fitted.values <- -fitted.values }else fitted.values <- pred weights_update <- .Fortran("update_wt_ccsvm", n=as.integer(n), weights=as.double(weights), y=as.double(y), f=as.double(fitted.values), cfun=as.integer(cfun), dfun=as.integer(dfun), s=as.double(s), eps=as.double(epsilon), delta=as.double(delta), weights_update=as.double(rep(0, n)), PACKAGE="mpath")$weights_update subset <- which(weights_update > 0) if(trace && kernel=="linear"){ beta <- t(RET$coefs) %*% as.matrix(xtmp[RET$index,]) #kernel="linear" only: convert dual coef to RET$coefs to primal coef #RET$b0 <- -RET$rho #intercept ###for some reason, depending on the sign of the first element of y, in wsvm, version 1.7-5 or svm version 1.7.3, beta and b0 are minus corresponding values coef_ and intercetp_ in Python scikit.learn.svm.SVC version xxx. Likewise, fitted values have inverse sign ### furthermore, it is questionable on the beta value from wsvm sometimes, but fitted values seem ok. if(dfun==6 && y[1]==-1) beta <- -beta los[k] <- .Fortran("loss2_ccsvm", n=as.integer(n), y=as.double(y), f=as.double(fitted.values), weights=as.double(weights), cfun=as.integer(cfun), dfun=as.integer(dfun), s=as.double(s), eps=as.double(epsilon), delta=as.double(delta), los=as.double(0.0), PACKAGE="mpath")$los los[k] <- los[k]+0.5/cost*sum(beta^2) ### only for kernel="linear" cat("\niteration", k, ", robust loss value", los[k], "\n") if(k > 1){ d1 <- abs((los[k]-los[k-1])/los[k-1]) if(los[k] > los[k-1]){ k <- iter } } }else{ if(k > 1){ d1 <- sum((fitold-fitted.values)^2)/sum(fitted.values^2) } fitold <- fitted.values } k <- k + 1 if(trace) cat("d1=", d1, ", k=", k, ", d1 > reltol && k <= iter: ", (d1 > reltol && k <= iter), "\n") } RET$cfun <- cfunold RET$s <- s RET$delta <- delta RET$weights_update <- weights_update RET }