./PaxHeaders.7814/EEL.R0000644000000000000000000000013112567312012012661 xustar000000000000000030 mtime=1440584714.746421615 29 atime=1440584630.65842735 30 ctime=1440584714.746421615 EEL.R0000644000175000001440000005000012567312012012322 0ustar00stsmalusers00000000000000./PaxHeaders.7808/EL-Owen.R0000644000000000000000000000013012567311666013502 xustar000000000000000029 mtime=1440584630.65942735 30 atime=1440584658.298425465 29 ctime=1440584630.65942735 EL-Owen.R0000644000175000001440000001631212567311666013151 0ustar00stsmalusers00000000000000# # The function elm() calculates the empirical likelihood # ratio for the mean. # # INPUTS: # x Matrix (or vector) of data # One row (or element) per observation # mu Value to be tested as mean for x observations # lam Optional starting value for Lagrange multiplier # maxit Optional maximum number of iterations # gradtol Optional tolerance for convergence test # svdtol Optional tolerance for detecting singularity # while solving equations # itertrace Optional flag to print results during iteration # # OUTPUTS: # logelr Log empirical likelihood # lambda Lagrange multiplier # grad Gradient of log likelihood # hess Hessian of log likelihood # wts Relative observation weights at the solution # nits Number of iterations used # If mu is in the interior of the convex hull of the # observations x, then wts should sum to n. If mu is outside # the convex hull then wts should sum to nearly zero, and logelr # will be a large negative number. It should be minus infinity, # but for inferential purposes a very large negative number is # essentially equivalent. If mu is on the boundary of the convex # hull then wts should sum to nearly k where k is the number of # observations within that face of the convex hull which contains mu. # # When mu is interior to the convex hull, it is typical for # the algorithm to converge quadratically to the solution, perhaps # after a few iterations of searching to get near the solution. # When mu is outside or near the boundary of the convex hull, then # the solution involves a lambda of infinite norm. The algorithm # tends to nearly double lambda at each iteration and the gradient # size then decreases roughly by half at each iteration. # # The goal in writing the algorithm was to have it "fail gracefully" # when mu is not inside the convex hull. The user can either leave # logelr "large and negative" or can replace it by minus infinity when # the weights do not sum to nearly n. # ## EL from Own el.owen <- function (x, mu, lam, maxit=25, gradtol=1e-7, svdtol=1e-9, itertrace=F, TINY=sqrt(.Machine$double.xmin)) { x <- as.matrix(x) n <- nrow(x) p <- ncol(x) mu <- as.vector(mu) if (length(mu) != p) stop("mu must have same dimension as observation vectors.") if (n <= p) stop("Need more observations than length(mu) in el.test().") z <- t(t(x) - mu) scale <- mean(abs(z)) + TINY z <- z/scale if (!missing(lam)) { lam <- as.vector(lam) lam <- lam * scale if (logelr(z, rep(0, p), lam) > 0) lam <- rep(0, p) } if (missing(lam)) lam <- rep(0, p) if (svdtol < TINY) svdtol <- TINY if (gradtol < TINY) gradtol <- TINY nwts <- c(3^-c(0:3), rep(0, 12)) gwts <- 2^(-c(0:(length(nwts) - 1))) gwts <- (gwts^2 - nwts^2)^0.5 gwts[12:16] <- gwts[12:16] * 10^-c(1:5) nits <- 0 gsize <- gradtol + 1 while (nits < maxit && gsize > gradtol) { arg <- 1 + z %*% lam wts1 <- as.vector(llogp(arg, 1/(n))) wts2 <- as.vector(-llogpp(arg, 1/(n)))^0.5 grad <- as.matrix(-z * wts1) grad <- as.vector(rowsum(grad, rep(1, nrow(grad)))) gsize <- mean(abs(grad)) hess <- z * wts2 svdh <- svd(hess) if (min(svdh$d) < max(svdh$d) * svdtol) svdh$d <- svdh$d + max(svdh$d) * svdtol nstep <- svdh$v %*% (t(svdh$u)/svdh$d) nstep <- as.vector(nstep %*% matrix(wts1/wts2, n, 1)) gstep <- -grad if (sum(nstep^2) < sum(gstep^2)) gstep <- gstep * (sum(nstep^2)^0.5/sum(gstep^2)^0.5) ologelr <- -sum(llog(arg, 1/(n))) ninner <- 0 for (i in 1:length(nwts)) { nlogelr <- logelr(z, rep(0, p), lam + nwts[i] * nstep + gwts[i] * gstep) if (nlogelr < ologelr) { lam <- lam + nwts[i] * nstep + gwts[i] * gstep ninner <- i break } } nits <- nits + 1 if (ninner == 0) nits <- maxit if (itertrace) print(c(lam, nlogelr, gsize, ninner)) } list(`-2LLR` = -2 * nlogelr, Pval = 1 - pchisq(-2 * nlogelr, df = p), lambda = lam/scale, grad = grad * scale, hess = t(hess) %*% hess * scale^2, wts = wts1, nits = nits) } logelr <- function (x, mu, lam) { x <- as.matrix(x) n <- nrow(x) p <- ncol(x) if (n <= p) stop("Need more observations than variables in logelr.") mu <- as.vector(mu) if (length(mu) != p) stop("Length of mean doesn't match number of variables in logelr.") z <- t(t(x) - mu) arg <- 1 + z %*% lam return(-sum(llog(arg, 1/n))) } llog <- function (z, eps) { ans <- z lo <- (z < eps) ans[lo] <- log(eps) - 1.5 + 2 * z[lo]/eps - 0.5 * (z[lo]/eps)^2 ans[!lo] <- log(z[!lo]) ans } llogp <- function (z, eps) { ans <- z lo <- (z < eps) ans[lo] <- 2/eps - z[lo]/eps^2 ans[!lo] <- 1/z[!lo] ans } llogpp <- function (z, eps) { ans <- z lo <- (z < eps) ans[lo] <- -1/eps^2 ans[!lo] <- -1/z[!lo]^2 ans } ## EL from Own, for estimating equation ## mx: m(theta,x) el.ee <- function (mx, lam, maxit=25, gradtol=1e-7, svdtol=1e-9, itertrace=F, TINY=sqrt(.Machine$double.xmin)) { mx <- as.matrix(mx) n <- nrow(mx) p <- ncol(mx) if (n <= p) stop("Need more observations than length(mu) in el.test().") z <- mx scale <- mean(abs(z)) + TINY z <- z/scale if (!missing(lam)) { lam <- as.vector(lam) lam <- lam * scale if (logelr.ee(z, lam) > 0) lam <- rep(0, p) } if (missing(lam)) lam <- rep(0, p) if (svdtol < TINY) svdtol <- TINY if (gradtol < TINY) gradtol <- TINY nwts <- c(3^-c(0:3), rep(0, 12)) gwts <- 2^(-c(0:(length(nwts) - 1))) gwts <- (gwts^2 - nwts^2)^0.5 gwts[12:16] <- gwts[12:16] * 10^-c(1:5) nits <- 0 gsize <- gradtol + 1 while (nits < maxit && gsize > gradtol) { arg <- 1 + z %*% lam wts1 <- as.vector(llogp(arg, 1/(n))) wts2 <- as.vector(-llogpp(arg, 1/(n)))^0.5 grad <- as.matrix(-z * wts1) grad <- as.vector(rowsum(grad, rep(1, nrow(grad)))) gsize <- mean(abs(grad)) hess <- z * wts2 svdh <- svd(hess) if (min(svdh$d) < max(svdh$d) * svdtol) svdh$d <- svdh$d + max(svdh$d) * svdtol nstep <- svdh$v %*% (t(svdh$u)/svdh$d) nstep <- as.vector(nstep %*% matrix(wts1/wts2, n, 1)) gstep <- -grad if (sum(nstep^2) < sum(gstep^2)) gstep <- gstep * (sum(nstep^2)^0.5/sum(gstep^2)^0.5) ologelr <- -sum(llog(arg, 1/(n))) ninner <- 0 for (i in 1:length(nwts)) { nlogelr <- logelr.ee(z, lam + nwts[i] * nstep + gwts[i] * gstep) if (nlogelr < ologelr) { lam <- lam + nwts[i] * nstep + gwts[i] * gstep ninner <- i break } } nits <- nits + 1 if (ninner == 0) nits <- maxit if (itertrace) print(c(lam, nlogelr, gsize, ninner)) } list(`-2LLR` = -2 * nlogelr, Pval = 1 - pchisq(-2 * nlogelr, df = p), lambda = lam/scale, grad = grad * scale, hess = t(hess) %*% hess * scale^2, wts = wts1, nits = nits) } logelr.ee <- function (mx, lam) { mx <- as.matrix(mx) n <- nrow(mx) p <- ncol(mx) if (n <= p) stop("Need more observations than variables in logelr.") z <- mx arg <- 1 + z %*% lam return(-sum(llog(arg, 1/n))) } ./PaxHeaders.7808/EL.R0000644000000000000000000000013012567311666012574 xustar000000000000000029 mtime=1440584630.65942735 30 atime=1440584700.372422596 29 ctime=1440584630.65942735 EL.R0000644000175000001440000001637312567311666012252 0ustar00stsmalusers00000000000000source('EL-Owen.R') ###################################################################### ## The first direvative of scad penalty, see Fan and Li (2001, Eq 2.7) ## Input: mu - the argument ###################################################################### scad.grad=function(mu, scadlambda, Dim=length(mu), scada=3.7){ mu=abs(mu) t2=ifelse((mu<=scadlambda),1,0) t3=(scada*scadlambda-mu) t3=ifelse(t3>=0,t3,0) t3=t3/((scada-1)*scadlambda) t3=ifelse(mu>scadlambda,t3,0) ss=scadlambda*(t2+t3) ss[mu==0]=0 ss } #################################################################### ## SCAD penalty function p(), se Fan and Li (2006, pp603) #################################################################### scad=function(mu, scadlambda, Dim=length(mu), scada=3.7){ mu=abs(mu) tmp=muscadlambda)&(mu<=scada*scadlambda) s[tmp]=-(mu[tmp]*mu[tmp]-2*scada*scadlambda*mu[tmp]+ scadlambda*scadlambda)/2/(scada-1) tmp=mu>=scada*scadlambda s[tmp]=(scada+1)*scadlambda*scadlambda/2 s } ################################################################ # Local quadratic approximation # Iterative procedure ################################################################ elhood.lqa=function(init.mu, x, scadlambda, maxiter=30){ n=nrow(x) p=ncol(x) ## determining zeros sd.mu=apply(x,2,sd) threshold=sd.mu*thres.lv ## criterion to determine zero # init value for lambda init.lambda=rep(0, p) #nzero=abs(init.mu)>threshold #init.mu[!nzero]=0 ## mu: the nonzero mu's, mu.all: all the mu's (zeo or nonzero) elhood=function(mu){ mu.all=rep(0,p) mu.all[!zero]=mu lambda=el.owen(x,mu.all,init.lambda)$lambda z=t(t(x)-mu.all) arg=1+z%*%lambda nlogelr=sum(llog(arg,1/n)) tmp=scad.grad(init.mu, scadlambda)/abs(init.mu) vf=nlogelr+sum(tmp*mu*mu)*n/2 vgm=tmp*mu*n-lambda[!zero]*sum(llogp(arg,1/n)) vf=nlogelr+sum(scad(mu,scadlambda))*n vgm=-lambda[!zero]*sum(llogp(arg,1/n))+n*scad.grad(mu,scadlambda)*sign(mu) attr(vf, "gradient")=vgm return(vf) } counter=1 init.mu.all=init.mu zero=abs(init.mu)p-1 | max(abs(dff))<1e-3) break #browser() #if(zz$code<4) # break #if(max(abs(zz$grad))<1e-3) # break #if(max(abs(dff))<1e-3) # break counter=counter+1 if (counter >= 1) { warning("iteration fails to converge") break } } #minimum value est=init.mu.all z=t(t(x)-est) arg=1+z%*%lambda val=2*sum(llog(arg,1/n)) #browser() list(estimate=est, val=val) } ########################################################## ## Evaluating the EL ratio for fixed components of mu ## loci is the location indicator vector, loci.v is the value given ############################################################ elhood.lqa.fix=function(init.mu, x, scadlambda, maxiter=30,loci,loci.v){ n=nrow(x) p=ncol(x) ## determining zeros sd.mu=apply(x,2,sd) threshold=sd.mu*thres.lv ## criterion to determine zero # init value for lambda init.lambda=rep(0, p) #nzero=abs(init.mu)>threshold #init.mu[!nzero]=0 ## mu: the nonzero mu's, mu.all: all the mu's (zeo or nonzero) elhood=function(mu){ mu.all=rep(0,p) mu.all[loci]=loci.v mu.all[!zero]=mu lambda=el.owen(x,mu.all,init.lambda)$lambda z=t(t(x)-mu.all) arg=1+z%*%lambda nlogelr=sum(llog(arg,1/n)) tmp=scad.grad(init.mu, scadlambda)/abs(init.mu) vf=nlogelr+sum(tmp*mu*mu)*n/2 vgm=tmp*mu*n-lambda[!zero]*sum(llogp(arg,1/n)) vf=nlogelr+sum(scad(mu,scadlambda))*n vgm=-lambda[!zero]*sum(llogp(arg,1/n))+n*scad.grad(mu,scadlambda)*sign(mu) attr(vf, "gradient")=vgm return(vf) } counter=1 init.mu.all=init.mu init.mu[loci]=0.0 zero=abs(init.mu)p-1 | max(abs(dff))<1e-3) break #browser() #if(zz$code<4) # break #if(max(abs(zz$grad))<1e-3) # break #if(max(abs(dff))<1e-3) # break counter=counter+1 if (counter >= 1) { warning("iteration fails to converge") break } } #minimum value est=init.mu.all z=t(t(x)-est) arg=1+z%*%lambda val=2*sum(llog(arg,1/n)) #browser() list(estimate=est, val=val) } ################################################################## # BIC for EL # See Wang, Li and Leng (JRSSB, 2009), Chen, J. and Chen, Z. (2008) # Input: mu - initial value # x - the data matrix # scadlambda - a vector of lambda in scad # Output: est - the fitted mu's, one row for each scad lambda # bic - bic values # bicc - bic in Wang, Li and Leng # ebic - the extended bic in Chen and Chen # est.bic - the estimated mu chosen by BIC # est.bicc - the estimated mu chosen by BIC in Wang, Li and Leng # est.ebic - the estimated mu chosen by EBIC ################################################################# bic.el=function(mu,x,scadlambda,N=nrow(x)){ nl=length(scadlambda) val=rep(0,nl) est=matrix(0,nl,ncol(x)) mu.init=mu for(i in 1:nl) { # cat(scadlambda[i],' ') fit=elhood.lqa(mu.init,x,scadlambda[i]) val[i]=fit$val est[i,]=fit$estimate mu.init=est[i,] } cat('\n') dof=apply(abs(est)>0,1,sum) bic=val+dof*log(N) bicc=val+dof*log(N)*max(1,log(log(N))) ebic=val+dof*(log(N)+2*log(ncol(x))) s.bic=sort(bic,index=TRUE) est.bic=est[s.bic$ix[1],] s.bicc=sort(bicc,index=TRUE) est.bicc=est[s.bicc$ix[1],] s.ebic=sort(ebic,index=TRUE) est.ebic=est[s.ebic$ix[1],] list(est=est,bic=bic,bicc=bicc,ebic=ebic,dof=dof, est.bic=est.bic,est.bicc=est.bicc,est.ebic=est.ebic,l1=s.bic$ix[1],l2=s.bicc$ix[1],l3=s.ebic$ix[i],vv=val) ### added location of the scad lambda } ./PaxHeaders.7814/EL-Owen.R0000644000000000000000000000013012567311666013477 xustar000000000000000029 mtime=1440584630.65942735 30 atime=1440584658.298425465 29 ctime=1440584630.65942735 EL-Owen.R0000644000175000001440000001631212567311666013151 0ustar00stsmalusers00000000000000# # The function elm() calculates the empirical likelihood # ratio for the mean. # # INPUTS: # x Matrix (or vector) of data # One row (or element) per observation # mu Value to be tested as mean for x observations # lam Optional starting value for Lagrange multiplier # maxit Optional maximum number of iterations # gradtol Optional tolerance for convergence test # svdtol Optional tolerance for detecting singularity # while solving equations # itertrace Optional flag to print results during iteration # # OUTPUTS: # logelr Log empirical likelihood # lambda Lagrange multiplier # grad Gradient of log likelihood # hess Hessian of log likelihood # wts Relative observation weights at the solution # nits Number of iterations used # If mu is in the interior of the convex hull of the # observations x, then wts should sum to n. If mu is outside # the convex hull then wts should sum to nearly zero, and logelr # will be a large negative number. It should be minus infinity, # but for inferential purposes a very large negative number is # essentially equivalent. If mu is on the boundary of the convex # hull then wts should sum to nearly k where k is the number of # observations within that face of the convex hull which contains mu. # # When mu is interior to the convex hull, it is typical for # the algorithm to converge quadratically to the solution, perhaps # after a few iterations of searching to get near the solution. # When mu is outside or near the boundary of the convex hull, then # the solution involves a lambda of infinite norm. The algorithm # tends to nearly double lambda at each iteration and the gradient # size then decreases roughly by half at each iteration. # # The goal in writing the algorithm was to have it "fail gracefully" # when mu is not inside the convex hull. The user can either leave # logelr "large and negative" or can replace it by minus infinity when # the weights do not sum to nearly n. # ## EL from Own el.owen <- function (x, mu, lam, maxit=25, gradtol=1e-7, svdtol=1e-9, itertrace=F, TINY=sqrt(.Machine$double.xmin)) { x <- as.matrix(x) n <- nrow(x) p <- ncol(x) mu <- as.vector(mu) if (length(mu) != p) stop("mu must have same dimension as observation vectors.") if (n <= p) stop("Need more observations than length(mu) in el.test().") z <- t(t(x) - mu) scale <- mean(abs(z)) + TINY z <- z/scale if (!missing(lam)) { lam <- as.vector(lam) lam <- lam * scale if (logelr(z, rep(0, p), lam) > 0) lam <- rep(0, p) } if (missing(lam)) lam <- rep(0, p) if (svdtol < TINY) svdtol <- TINY if (gradtol < TINY) gradtol <- TINY nwts <- c(3^-c(0:3), rep(0, 12)) gwts <- 2^(-c(0:(length(nwts) - 1))) gwts <- (gwts^2 - nwts^2)^0.5 gwts[12:16] <- gwts[12:16] * 10^-c(1:5) nits <- 0 gsize <- gradtol + 1 while (nits < maxit && gsize > gradtol) { arg <- 1 + z %*% lam wts1 <- as.vector(llogp(arg, 1/(n))) wts2 <- as.vector(-llogpp(arg, 1/(n)))^0.5 grad <- as.matrix(-z * wts1) grad <- as.vector(rowsum(grad, rep(1, nrow(grad)))) gsize <- mean(abs(grad)) hess <- z * wts2 svdh <- svd(hess) if (min(svdh$d) < max(svdh$d) * svdtol) svdh$d <- svdh$d + max(svdh$d) * svdtol nstep <- svdh$v %*% (t(svdh$u)/svdh$d) nstep <- as.vector(nstep %*% matrix(wts1/wts2, n, 1)) gstep <- -grad if (sum(nstep^2) < sum(gstep^2)) gstep <- gstep * (sum(nstep^2)^0.5/sum(gstep^2)^0.5) ologelr <- -sum(llog(arg, 1/(n))) ninner <- 0 for (i in 1:length(nwts)) { nlogelr <- logelr(z, rep(0, p), lam + nwts[i] * nstep + gwts[i] * gstep) if (nlogelr < ologelr) { lam <- lam + nwts[i] * nstep + gwts[i] * gstep ninner <- i break } } nits <- nits + 1 if (ninner == 0) nits <- maxit if (itertrace) print(c(lam, nlogelr, gsize, ninner)) } list(`-2LLR` = -2 * nlogelr, Pval = 1 - pchisq(-2 * nlogelr, df = p), lambda = lam/scale, grad = grad * scale, hess = t(hess) %*% hess * scale^2, wts = wts1, nits = nits) } logelr <- function (x, mu, lam) { x <- as.matrix(x) n <- nrow(x) p <- ncol(x) if (n <= p) stop("Need more observations than variables in logelr.") mu <- as.vector(mu) if (length(mu) != p) stop("Length of mean doesn't match number of variables in logelr.") z <- t(t(x) - mu) arg <- 1 + z %*% lam return(-sum(llog(arg, 1/n))) } llog <- function (z, eps) { ans <- z lo <- (z < eps) ans[lo] <- log(eps) - 1.5 + 2 * z[lo]/eps - 0.5 * (z[lo]/eps)^2 ans[!lo] <- log(z[!lo]) ans } llogp <- function (z, eps) { ans <- z lo <- (z < eps) ans[lo] <- 2/eps - z[lo]/eps^2 ans[!lo] <- 1/z[!lo] ans } llogpp <- function (z, eps) { ans <- z lo <- (z < eps) ans[lo] <- -1/eps^2 ans[!lo] <- -1/z[!lo]^2 ans } ## EL from Own, for estimating equation ## mx: m(theta,x) el.ee <- function (mx, lam, maxit=25, gradtol=1e-7, svdtol=1e-9, itertrace=F, TINY=sqrt(.Machine$double.xmin)) { mx <- as.matrix(mx) n <- nrow(mx) p <- ncol(mx) if (n <= p) stop("Need more observations than length(mu) in el.test().") z <- mx scale <- mean(abs(z)) + TINY z <- z/scale if (!missing(lam)) { lam <- as.vector(lam) lam <- lam * scale if (logelr.ee(z, lam) > 0) lam <- rep(0, p) } if (missing(lam)) lam <- rep(0, p) if (svdtol < TINY) svdtol <- TINY if (gradtol < TINY) gradtol <- TINY nwts <- c(3^-c(0:3), rep(0, 12)) gwts <- 2^(-c(0:(length(nwts) - 1))) gwts <- (gwts^2 - nwts^2)^0.5 gwts[12:16] <- gwts[12:16] * 10^-c(1:5) nits <- 0 gsize <- gradtol + 1 while (nits < maxit && gsize > gradtol) { arg <- 1 + z %*% lam wts1 <- as.vector(llogp(arg, 1/(n))) wts2 <- as.vector(-llogpp(arg, 1/(n)))^0.5 grad <- as.matrix(-z * wts1) grad <- as.vector(rowsum(grad, rep(1, nrow(grad)))) gsize <- mean(abs(grad)) hess <- z * wts2 svdh <- svd(hess) if (min(svdh$d) < max(svdh$d) * svdtol) svdh$d <- svdh$d + max(svdh$d) * svdtol nstep <- svdh$v %*% (t(svdh$u)/svdh$d) nstep <- as.vector(nstep %*% matrix(wts1/wts2, n, 1)) gstep <- -grad if (sum(nstep^2) < sum(gstep^2)) gstep <- gstep * (sum(nstep^2)^0.5/sum(gstep^2)^0.5) ologelr <- -sum(llog(arg, 1/(n))) ninner <- 0 for (i in 1:length(nwts)) { nlogelr <- logelr.ee(z, lam + nwts[i] * nstep + gwts[i] * gstep) if (nlogelr < ologelr) { lam <- lam + nwts[i] * nstep + gwts[i] * gstep ninner <- i break } } nits <- nits + 1 if (ninner == 0) nits <- maxit if (itertrace) print(c(lam, nlogelr, gsize, ninner)) } list(`-2LLR` = -2 * nlogelr, Pval = 1 - pchisq(-2 * nlogelr, df = p), lambda = lam/scale, grad = grad * scale, hess = t(hess) %*% hess * scale^2, wts = wts1, nits = nits) } logelr.ee <- function (mx, lam) { mx <- as.matrix(mx) n <- nrow(mx) p <- ncol(mx) if (n <= p) stop("Need more observations than variables in logelr.") z <- mx arg <- 1 + z %*% lam return(-sum(llog(arg, 1/n))) } ./PaxHeaders.7814/EL.R0000644000000000000000000000013012567311666012571 xustar000000000000000029 mtime=1440584630.65942735 30 atime=1440584700.372422596 29 ctime=1440584630.65942735 EL.R0000644000175000001440000001637312567311666012252 0ustar00stsmalusers00000000000000source('EL-Owen.R') ###################################################################### ## The first direvative of scad penalty, see Fan and Li (2001, Eq 2.7) ## Input: mu - the argument ###################################################################### scad.grad=function(mu, scadlambda, Dim=length(mu), scada=3.7){ mu=abs(mu) t2=ifelse((mu<=scadlambda),1,0) t3=(scada*scadlambda-mu) t3=ifelse(t3>=0,t3,0) t3=t3/((scada-1)*scadlambda) t3=ifelse(mu>scadlambda,t3,0) ss=scadlambda*(t2+t3) ss[mu==0]=0 ss } #################################################################### ## SCAD penalty function p(), se Fan and Li (2006, pp603) #################################################################### scad=function(mu, scadlambda, Dim=length(mu), scada=3.7){ mu=abs(mu) tmp=muscadlambda)&(mu<=scada*scadlambda) s[tmp]=-(mu[tmp]*mu[tmp]-2*scada*scadlambda*mu[tmp]+ scadlambda*scadlambda)/2/(scada-1) tmp=mu>=scada*scadlambda s[tmp]=(scada+1)*scadlambda*scadlambda/2 s } ################################################################ # Local quadratic approximation # Iterative procedure ################################################################ elhood.lqa=function(init.mu, x, scadlambda, maxiter=30){ n=nrow(x) p=ncol(x) ## determining zeros sd.mu=apply(x,2,sd) threshold=sd.mu*thres.lv ## criterion to determine zero # init value for lambda init.lambda=rep(0, p) #nzero=abs(init.mu)>threshold #init.mu[!nzero]=0 ## mu: the nonzero mu's, mu.all: all the mu's (zeo or nonzero) elhood=function(mu){ mu.all=rep(0,p) mu.all[!zero]=mu lambda=el.owen(x,mu.all,init.lambda)$lambda z=t(t(x)-mu.all) arg=1+z%*%lambda nlogelr=sum(llog(arg,1/n)) tmp=scad.grad(init.mu, scadlambda)/abs(init.mu) vf=nlogelr+sum(tmp*mu*mu)*n/2 vgm=tmp*mu*n-lambda[!zero]*sum(llogp(arg,1/n)) vf=nlogelr+sum(scad(mu,scadlambda))*n vgm=-lambda[!zero]*sum(llogp(arg,1/n))+n*scad.grad(mu,scadlambda)*sign(mu) attr(vf, "gradient")=vgm return(vf) } counter=1 init.mu.all=init.mu zero=abs(init.mu)p-1 | max(abs(dff))<1e-3) break #browser() #if(zz$code<4) # break #if(max(abs(zz$grad))<1e-3) # break #if(max(abs(dff))<1e-3) # break counter=counter+1 if (counter >= 1) { warning("iteration fails to converge") break } } #minimum value est=init.mu.all z=t(t(x)-est) arg=1+z%*%lambda val=2*sum(llog(arg,1/n)) #browser() list(estimate=est, val=val) } ########################################################## ## Evaluating the EL ratio for fixed components of mu ## loci is the location indicator vector, loci.v is the value given ############################################################ elhood.lqa.fix=function(init.mu, x, scadlambda, maxiter=30,loci,loci.v){ n=nrow(x) p=ncol(x) ## determining zeros sd.mu=apply(x,2,sd) threshold=sd.mu*thres.lv ## criterion to determine zero # init value for lambda init.lambda=rep(0, p) #nzero=abs(init.mu)>threshold #init.mu[!nzero]=0 ## mu: the nonzero mu's, mu.all: all the mu's (zeo or nonzero) elhood=function(mu){ mu.all=rep(0,p) mu.all[loci]=loci.v mu.all[!zero]=mu lambda=el.owen(x,mu.all,init.lambda)$lambda z=t(t(x)-mu.all) arg=1+z%*%lambda nlogelr=sum(llog(arg,1/n)) tmp=scad.grad(init.mu, scadlambda)/abs(init.mu) vf=nlogelr+sum(tmp*mu*mu)*n/2 vgm=tmp*mu*n-lambda[!zero]*sum(llogp(arg,1/n)) vf=nlogelr+sum(scad(mu,scadlambda))*n vgm=-lambda[!zero]*sum(llogp(arg,1/n))+n*scad.grad(mu,scadlambda)*sign(mu) attr(vf, "gradient")=vgm return(vf) } counter=1 init.mu.all=init.mu init.mu[loci]=0.0 zero=abs(init.mu)p-1 | max(abs(dff))<1e-3) break #browser() #if(zz$code<4) # break #if(max(abs(zz$grad))<1e-3) # break #if(max(abs(dff))<1e-3) # break counter=counter+1 if (counter >= 1) { warning("iteration fails to converge") break } } #minimum value est=init.mu.all z=t(t(x)-est) arg=1+z%*%lambda val=2*sum(llog(arg,1/n)) #browser() list(estimate=est, val=val) } ################################################################## # BIC for EL # See Wang, Li and Leng (JRSSB, 2009), Chen, J. and Chen, Z. (2008) # Input: mu - initial value # x - the data matrix # scadlambda - a vector of lambda in scad # Output: est - the fitted mu's, one row for each scad lambda # bic - bic values # bicc - bic in Wang, Li and Leng # ebic - the extended bic in Chen and Chen # est.bic - the estimated mu chosen by BIC # est.bicc - the estimated mu chosen by BIC in Wang, Li and Leng # est.ebic - the estimated mu chosen by EBIC ################################################################# bic.el=function(mu,x,scadlambda,N=nrow(x)){ nl=length(scadlambda) val=rep(0,nl) est=matrix(0,nl,ncol(x)) mu.init=mu for(i in 1:nl) { # cat(scadlambda[i],' ') fit=elhood.lqa(mu.init,x,scadlambda[i]) val[i]=fit$val est[i,]=fit$estimate mu.init=est[i,] } cat('\n') dof=apply(abs(est)>0,1,sum) bic=val+dof*log(N) bicc=val+dof*log(N)*max(1,log(log(N))) ebic=val+dof*(log(N)+2*log(ncol(x))) s.bic=sort(bic,index=TRUE) est.bic=est[s.bic$ix[1],] s.bicc=sort(bicc,index=TRUE) est.bicc=est[s.bicc$ix[1],] s.ebic=sort(ebic,index=TRUE) est.ebic=est[s.ebic$ix[1],] list(est=est,bic=bic,bicc=bicc,ebic=ebic,dof=dof, est.bic=est.bic,est.bicc=est.bicc,est.ebic=est.ebic,l1=s.bic$ix[1],l2=s.bicc$ix[1],l3=s.ebic$ix[i],vv=val) ### added location of the scad lambda }