fast.s <- function(x, y, int, N, k=2, best.r=5, b=.5, cc=1.56, seed) { # # int = 1 -> add a column of ones to x # N = cant de sub-samples # k = number of refining iterations in ea. # subsample # k = 0 means "raw-subsampling" # b = right hand side of the equation # cc = corresponding tuning constant # best.r = number of "best betas" to remember # from the subsamples. These will be later # iterated until convergence # the objective function, we solve loss.S(u,s,cc)=b for "s" loss.S <- function(u,s,cc) mean(rho(u/s,cc) ) re.s <- function(x,y,initial.beta,initial.scale,k,conv,b,cc) { # does "k" IRWLS refining steps from "initial.beta" # # if "initial.scale" is present, it's used, o/w the MAD is used # k = number of refining steps # conv = 0 means "do k steps and don't check for convergence" # conv = 1 means "stop when convergence is detected, or the # maximum number of iterations is achieved" # b and cc = tuning constants of the equation # f.w <- function(u, cc) { # weight function = psi(u)/u tmp <- (1 - (u/cc)^2)^2 tmp[ abs(u/cc) > 1 ] <- 0 return(tmp) } n <- dim(x)[1] p <- dim(x)[2] res <- y - x %*% initial.beta if( missing( initial.scale ) ) initial.scale <- scale <- median(abs(res))/.6745 else scale <- initial.scale if( conv == 1) k <- 50 # # if conv == 1 then set the max no. of iterations to 50 # magic number alert!!! beta <- initial.beta lower.bound <- median(abs(res))/cc for(i in 1:k) { # do one step of the iterations to solve for the scale scale.super.old <- scale #lower.bound <- median(abs(res))/1.56 scale <- sqrt( scale^2 * mean( rho( res / scale, cc ) ) / b ) # now do one step of IRWLS with the "improved scale" weights <- f.w( res/scale, cc ) W <- matrix(weights, n, p) xw <- x * sqrt(W) yw <- y * sqrt(weights) beta.1 <- our.solve( t(xw) %*% xw ,t(xw) %*% yw ) if(any(is.na(beta.1))) { beta.1 <- initial.beta scale <- initial.scale break } if( (conv==1) ) { # check for convergence if( norm( beta - beta.1 ) / norm(beta) < 1e-20 ) break # magic number alert!!! } res <- y - x %*% beta.1 beta <- beta.1 } res <- y - x %*% beta # get the residuals from the last beta return(list(beta.rw = beta.1, scale.rw = scale)) } scale1 <- function(u, b, cc, initial.sc=median(abs(u))/.6745) { # find the scale, full iterations max.it <- 200 # magic number alert #sc <- median(abs(u))/.6745 sc <- initial.sc i <- 0 eps <- 1e-20 # magic number alert err <- 1 while( ( (i <- i+1) < max.it ) && (err > eps) ) { sc2 <- sqrt( sc^2 * mean( rho( u / sc, cc ) ) / b ) err <- abs(sc2/sc - 1) sc <- sc2 } return(sc) } n <- dim(x)[1] p <- dim(x)[2] if( int == 1) { x <- cbind(rep(1,n), x) p <- p + 1 } if(!missing(seed)) set.seed(seed) best.betas <- matrix(0, best.r, p) best.scales <- rep(1e20, best.r) s.worst <- 1e20 n.ref <- 1 for(i in 1:N) { # # get a subsample # singular <- T while(singular==T) { indices <- sample(n,p) xs <- x[indices,] ys <- y[indices] beta <- our.solve(xs,ys) singular <- any(is.na(beta)) } if(k>0) { # do the refining tmp <- re.s(x=x,y=y,initial.beta=beta,k=k,conv=0,b=b,cc=cc) beta.rw <- tmp\$beta.rw scale.rw <- tmp\$scale.rw res.rw <- y - x %*% beta.rw } else { #k = 0 means "no refining" beta.rw <- beta res.rw <- y - x %*% beta.rw scale.rw <- median(abs(res.rw))/.6745 } if( i > 1 ) { # if this isn't the first iteration.... scale.test <- loss.S(res.rw,s.worst,cc) if( scale.test < b ) { s.best <- scale1(res.rw,b,cc,scale.rw) ind <- order(best.scales)[best.r] best.scales[ind] <- s.best best.betas[ind,] <- beta.rw s.worst <- max(best.scales) } } else { # if this is the first iteration, then this is the best beta... best.scales[best.r] <- scale1(res.rw,b,cc,scale.rw) best.betas[best.r,] <- beta.rw } } # do the complete refining step until convergence (conv=1) starting # from the best subsampling candidate (possibly refined) super.best.scale <- 1e20 # magic number alert for(i in best.r:1) { tmp <- re.s(x=x,y=y,initial.beta=best.betas[i,], initial.scale=best.scales[i],k=0,conv=1,b=b,cc=cc) if(tmp\$scale.rw < super.best.scale) { super.best.scale <- tmp\$scale.rw super.best.beta <- tmp\$beta.rw } } return(list(as.vector(super.best.beta), super.best.scale)) } our.solve <- function(a,b) { a <- qr(a) da <- dim(a\$qr) if(a\$rank < (p <- da[2])) return(NA) else qr.coef(a, b) } rho <- function(u, cc=1.56) { w <- abs(u)<=cc v <- (u^2/(2)*(1-(u^2/(cc^2))+(u^4/(3*cc^4))))*w +(1-w)*(cc^2/6) v <- v*6/cc^2 return(v) } norm <- function(x) sqrt( sum( x^2 ) )