#	R version June 16, 2004

# 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 <- cbind(X)}; 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 }
if(br>1) {br <<- br-1} else {browser()}
	ei <- y - a - X%*%b
	w^2/(2*sqrt(pi)*sig) - 2*w*mean(dnorm(ei,0,sig))
    }

    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)) }

    ans <- nlm(reg.w.mat.crit,x0,iterlim=100,print.level=pl)

    a <- ans$est[1]; b <- ans$est[2:(p+1)]; sig <- exp(ans$est[p+2])
    if(w.opt) {w <- exp(ans$est[p+3])} else {w <- win}

    invisible(list(a=a,b=b,sig=sig,w=w,res=c(y-a-X%*%b)))
}