# Copyright David W. Scott February 5, 2004 # Partial L2E multivariate linear regression # Purpose: fit model y = a + t(b)*x + e by L2E # the error distribution can be either N(0,sig^2) or w N(0,sig^2) # Inputs: # x - n x p input data matrix (p=1 vector OK) # y - n responses # xin - input guesses for (a,b,sig) length 1+p+1 # defaults = c(mean(y),0,..,0,sd(y)) # w.opt - False (default e ~ N(0,sig^2) or True (e ~ w N(0,sig^2)) # win - fixed value for w or initial guess for w (see w.opt) # pl - print level for nlmin (1=default 0=none) # Output: list containing estimated coefficients # \$a - intercept # \$b - vector of estimated beta's # \$sig - estimated standard deviation of residuals # \$w - estimated weight (note: w can be greater than 1) # \$res - residual vector reg.w.mat <- function(x,y,xin,w.opt=F,win=1,pl=1) { if(!is.matrix(x)) {x <- matrix(x,length(x),1)}; p <- ncol(x) if(length(y)!=nrow(x)) stop("x and y differ on n") reg.w.mat.crit <- function(x) { a <- x[1]; b <- x[2:(p+1)]; sig <- exp(x[p+2]) if(w.opt) { w <- exp(x[p+3]) } else { w <- win } ei <- y - a - X%*%b w^2/(2*sqrt(pi)*sig) - 2*w*mean(dnorm(ei,0,sig)) } assign("p",p,f=1); assign("w.opt",w.opt,f=1) assign("X",x,f=1); assign("y",y,f=1) if(missing(xin)) { xin <- c(mean(y),rep(0,p),sqrt(var(y))) } if(length(xin)!=p+2) stop("xin error") x0 <- c(xin[1:(p+1)],log(xin[p+2])) if(w.opt) { x0 _ c(x0,log(win)) } else { assign("win",win,f=1) } ans <- nlmin(reg.w.mat.crit,x0,max.iter=100,max.fcal=150,print.level=pl) a <- ans\$x[1]; b <- ans\$x[2:(p+1)]; sig <- exp(ans\$x[p+2]) if(w.opt) {w <- exp(ans\$x[p+3])} else {w<-win} invisible(list(a=a,b=b,sig=sig,w=w,res=c(y-a-x%*%b))) }