R/loss_rob.R

In nlr: Nonlinear Regression Modelling using Robust Methods

#************************************************************************
#**   robust loss unction                                              **
#**     start: argument must be a list containing "sigma".             **
#**     argument ... is the arguments for 'robfunc' extra argument     **
#**     for 'robfunc' the robust rho function.                         **
#**     Fault: = 9,  sigma is not a part of start list.                **
#************************************************************************
robloss <- function(formula,data,start,robfunc,control=nlr.control(robscale=T),...){
	Fault2 <- Fault()
	data <- as.list(data)
	start <- as.list(start)
	.tmp <-  (is.null(start$sigma))
	if(.tmp) return(list(Fault=Fault(FL=T,FN=8,FF="robloss")))
	datalist <- c(start,data)
	fmod <- eval(formula,datalist)						#  f(t) n*1									#
	if(is.Fault(fmod)) return(fmod)
	prd <- as.numeric(fmod$predictor)					#  f(t) values								#
	n <- length(prd)
	p <- formula$p
	gprd <- attr(fmod$predictor,"gradient")			#  f'(t)	Gradient	n*p						#
	hprd <- attr(fmod$predictor,"hessian")			#  f"(t)	Hessian	n*p*p					#
	rsd <- as.numeric(fmod$response) - prd			#  ri		Residals	n*1						#
  ris <- rsd / start$sigma 
	if(control$robscale)  robvalue <- robfunc$fnc(ris/robfunc$arguments$k1, ...)
	else robvalue <- robfunc$fnc(ris, ...)			#  rho(ri/sigma)	n*1							#
	rho <- as.numeric(robvalue)
	htheta <- sum(rho)
	grho <- attr(robvalue,"gradient")					#  V=rho'(ri/tht)		n*1						#
	hrho <- attr(robvalue,"hessian")					#  rho"(ri/tht)			n*1					#
	gres <- -1 * gprd										#  J=r'i=-grad(f)		n*p						#
	hres <- -hprd											#  r"i					n*p*p					#
	temp <- t(gres)
	gradh <- temp %*% grho								#  J'V, gradient of rho loss p*p			#
	gradh <- gradh / start$sigma						#  g(h) = 1/sgm * J'* V						#
	if(control$robscale) gradh <- gradh / robfunc$arguments$k1 
	dt <- diag(hrho)										#  D(thta)=diag[  rho"(ri/sgm)  ]   n*1	#
	hessh.p1 <- t(gres) %*% dt							#  J'D											#
	hessh.p1 <- hessh.p1 %*% gres						#  J'DJ										#
	hessh.p1 <- hessh.p1 / ( start$sigma ^ 2 )		#  hessian part 1 = 1/sigma^2  *  J'*D*J	#
	if(control$robscale) hessh.p1 <- hessh.p1 / robfunc$arguments$k1 ^ 2
	atheta <- matrix(rep(0,p*p), nrow=p)
	for(i in 1:n){
		a <- grho[i] * hres[i,,]
		atheta <- a + atheta
	}
	atheta <- atheta / start$sigma
	if(control$robscale) atheta <- atheta / robfunc$arguments$k1
	hessh.p2 <- atheta
	hessh <- hessh.p1 + hessh.p2
	dtilda <- diag( start$sigma * grho / rsd )		#  Dtilda(hteta) = rho' / (ri/sigma)		#
	attr(htheta,"gradient") <- gradh
	attr(htheta,"hessian") <- hessh
	attr(rsd,"gradient") <- gres
	attr(rsd,"hessian") <- hres
	result <- list(htheta=htheta,rho=robvalue,ri=rsd,hessh.p1=hessh.p1,hessh.p2=hessh.p2,dtilda=dtilda,fmod=fmod,Fault=Fault2)



	#****************************************************
	#*****   output:                                ***** 
	#*****  htheta=sum(rho)  fixed number           *****
	#*****    contains "gradient" & "hessian" attr  *****
	#*****  rho(ri/sigma)       n*1 vector          *****
	#*****    contains "gradient" & "hessian" attr  *****
	#*****  ri               residuals              *****
	#*****    contains "gradient" & "hessian" attr  *****
	#*****  h(theta)=he	ss.... p1+p2 hessian of (h)*****
	#*****  hess.p1 = p1                            *****
	#*****  hess.p2=  p2                            *****
	#*****  dtilda, new D(thilda) function         *****
	#*****  fmod:  computed function contains       *****
	#*****  response and or its gradient and hessian*****
	#*****  predictor and or its gradient & hessian *****
	#****************************************************

	return(result)
#	return(rho)
}