Nothing
#******************************************************************************
#* +------------------------------------------------------------------------+ *
#* | Function 'nlmest', M-estimate of a nonlinear function. | *
#* | Levenberg Marquardt Method. | *
#* | Note: becarefull to using this function when there is not outlier, it | *
#* | may not work witout outlier, in this case better to use nlmest | *
#* | the problem is in part of two (p2) in hessian its big here. | *
#* | argumnts: | *
#* | formula: 'nl.form' object, the function mode. | *
#* | data: data, contains dependents and independents, | *
#* | data.frame, list or named matrix it can be. | *
#* | start: starting values, it must contains 'sigma', selstart | *
#* | for nl.form object is not created yet, take cre of it. | *
#* | robfunc: obust function, it must be functionformat untill know, | *
#* | not a nl.form, in feature it should be modfied. | *
#* | ...: can be entries for robust loss function parameters. | *
#* | | *
#* +------------------------------------------------------------------------+ *
#******************************************************************************
nlmest.LM<-function(formula,data,start=getInitial(formula,data),robfunc,
control=nlr.control(tolerance=0.001, minlanda=1 / 2 ^ 10, maxiter=25 * length(start),robscale=T),vm=NULL,rm=eiginv(t(chol(vm))),...)
{
tolerance <- control$tolerance
maxiter <- control$maxiter
minlanda <- control$minlanda
trace <- control$trace
Fault2 <- Fault() ###### no error, automatically will be created by the Fault object
.tmp <- (is.null(start$sigma))
if (.tmp){
dt1 <- c(data,start)
names(dt1) <- c(names(data),names(start))
ht <- eval(formula,dt1)
if (is.Fault(ht)) return(ht)
dv <- as.numeric(ht$predictor) - as.numeric(ht$response)
if(! is.null(vm)) dv <- rm %*% dv
start[["sigma"]] <- mad(dv) # mscale(dv)
}
th <- start ## with sigma
theta <- unlist(start[names(start)!="sigma"]) ## without sigma
theta1 <- theta ## without sigma
datalist<-as.list(data)
.parameters <- names(start)
p <- length(theta)
n <- length(data[[1]])
datalist[.parameters] <- start[.parameters] #*****datalist contains both parameter vectors and data values
eol <- F
iterate <- 0
iterhist <- NULL
if(is.null(vm)) ht <- robloss(formula,data,th,robfunc,control=control,...) ## th: with sigma
else ht <- robloss.gn(formula,data,th,robfunc,rm,control=control,...)
if(any(is.missing(ht$ri))) {
return(nl.fitt.rob(Fault=Fault(FN=14,FF="nlmest.LM")))}
tmp <- as.numeric(ht$ri)
if(control$robscale) sigma <- mscale(tmp) #nl.mscale(tmp,robfunc,...)
else sigma <- sum(abs(tmp)) / n
if(is.Fault(sigma)) return(sigma)
#lines(data$xr,ht$fmod$predictor)
if(ht$Fault$FL) return(nl.fitt.rob(Fault=ht$Fault))
yresp <- as.numeric(ht$fmod$response)
ybar <- mean(yresp)
lambda <- 1
#///*************************************************
while (!eol) #///********* Start Iteration ****************
{
iterate <- iterate + 1
grh <- attr(ht$htheta,"gradient") ## g(theta) grad
hsh <- attr(ht$htheta,"hessian") ## H(theta) hess
ilev <-0
###################################################################
### LM part
###################################################################
heig <-eigen(hsh,symmetric=T) ## eig(f+land I) = (eig(F) +lnd)
lambda <- max(diag(hsh))*2 #heig$values)) * 2 ## add the smallest minus eig to
## to make all possitive
repeat{
ilev<- ilev + 1
I <- diag(dim(hsh) [2])
zeta <- hsh + lambda * diag(abs(diag(hsh)))
zetainv <- eiginv(zeta,stp=F,symmetric=T)
if(is.Fault(zetainv)) return(nl.fitt.rob(Fault=zetainv))
delta2 <- zetainv %*% grh
theta2 <- theta1 - delta2 #### without sigma
th2 <- as.list(theta2)
names(th2) <- names(formula$par)
th2["sigma"] <- sigma
if(is.null(vm)) ht2 <- robloss(formula,data,th2,robfunc,control=control,...)
else ht2 <- robloss.gn(formula,data,th2,robfunc,rm,control=control,...)
if(is.Fault(ht2)) return(nl.fitt.rob(Fault=ht2$Fault))
cnv <- sum((theta1-theta2)^2)
if(as.numeric(ht2$htheta) <= as.numeric(ht$htheta) || cnv < tolerance/10^10)
{
theta1 <- theta2
ht <- ht2
lamnbda <- lambda/10
break
}
else {
lambda <- lambda * 10
if( lambda > 1e12) return(nl.fitt.rob(Fault=Fault(FN=14,FF="nlmest.LM")))
}
}
g2 <- attr(ht$rho,"gradient") ## V(heta) = [rho.(r1/sg) ..... rho.(rn/sg)]T
tmp <- as.numeric(ht$ri)
if(control$robscale) sigma <- mscale(tmp) # nl.mscale(tmp,robfunc,...)
else sigma <- sum(abs(tmp)) / n
if(is.Fault(sigma)) return(sigma)
.expr1 <- t(attr(ht$ri,"gradient")) %*% attr(ht$ri,"gradient") ## J' J p*p
.expr2 <- eiginv(.expr1,symmetric=T,stp=F)
if(is.Fault(.expr2)) return(nl.fitt.rob(Fault=.expr2)) ## Note: eiginv only shoud return back Fault
.expr3 <- attr(ht$ri,"gradient") %*% .expr2 ## J (J' J)^-1 n*p
.expr4 <- .expr3 %*% t(attr(ht$ri,"gradient")) ## J (J' J)^-1 J' n*n
.expr5 <- .expr4 %*% g2 ## VT = J (J' J)^-1 J' V n*1
angle <- t(g2) %*% .expr5 ## V' * VT 1*1
angle <- angle / sqrt( sum(g2)^2 * sum(.expr5^2) )
th <-as.list(theta1) ## without sigma
names(th) <- names(formula$par)
th["sigma"] <- sigma ## renew with sigma
iterhist <- rbind(iterhist,c(iteration = iterate,objfnc = as.numeric(ht$htheta),
unlist(th),converge = angle ,ilev=ilev))
if(angle < tolerance) eol <- T
else if(iterate > maxiter) {
eol <- T
Fault2 <- Fault(FN=1,FF = "nlmest")
}
else {
if(is.null(vm)) ht <- robloss(formula,data,th,robfunc,control=control,...) ## h(theta) = sum( rho(ri/sigma) )
else ht <- robloss.gn(formula,data,th,robfunc,rm,control=control,...)
# lines(data$xr,ht$fmod$predictor)
if(ht$Fault$FL) return(nl.fitt.rob(Fault=ht$Fault))
}
} #\\\********* End Iteration ****************
#\\\*************************************************
if(! is.null(vm)){
rinv <- eiginv(rm)
yresp <- rinv %*% as.numeric(ht$fmod$response)
ypred <- rinv %*% as.numeric(ht$fmod$predictor)
}
else{
yresp <- as.numeric(ht$fmod$response)
ypred <- as.numeric(ht$fmod$predictor)
}
ybar <- mean(yresp)
nlrho <- 1 - sum( ( yresp - ypred ) ^ 2 ) /
sum( (ypred - ybar ) ^ 2 )
curv1 <- curvature(gradient = attr(ht$fmod$predictor,"gradient"),
hessian = attr(ht$fmod$predictor,"hessian"),
sigma = th[["sigma"]])
# curv2 <- curvature(gradient = attr(ht$fmod$predictor,"gradient"),
# hessian = attr(ht$fmod$predictor,"hessian"),
# sigma = th[["sigma"]])
if(is.null(vm))
result <- nl.fitt.rob(parameters = th,
scale = sigma,
correlation = nlrho,
form = formula,
response = ht$fmod$response,
predictor = ht$fmod$predictor,
curvature = curv1,
history = iterhist,
method = fittmethod(methodID= 7,
methodBR= switch(control$robscale,8,9), ### rob or nopnrob variance
detailBR= "Lev-Marq",
subroutine= "nlmest.LM"),
data = as.list(data),
sourcefnc = match.call(),
Fault = Fault2,
htheta = ht$htheta,
rho = ht$rho,
ri = ht$ri,
curvrob = NULL,
robform = robfunc
)
else
result <- nl.fitt.rgn(parameters = th,
scale = sigma,
correlation = nlrho,
form = formula,
response = transforminv(ht$fmod$response,rm),
predictor = transforminv(ht$fmod$predictor,rm),
curvature = curv1,
history = iterhist,
method = fittmethod(methodID= 7,
methodBR= switch(control$robscale,8,9), ### rob or nopnrob variance
detailBR= "Lev-Marq",
subroutine= "nlmest.LM"),
data = as.list(data),
sourcefnc = match.call(),
Fault = Fault2,
htheta = ht$htheta,
rho = ht$rho,
ri = ht$ri,
curvrob = NULL,
robform = robfunc,
vm = vm,
rm= rm,
gresponse = ht$fmod$response,
gpredictor = ht$fmod$predictor
)
return(result)
}
#******************************************************************************
#* +------------------------------------------------------------------------+ *
#* | End of 'nlmest.NLM' | *
#* | | *
#* | Hossein Riazoshams, UPM, INSPEM | *
#* | | *
#* | Writen Sep 2008, revised Mar 2009 | *
#* | | *
#* +------------------------------------------------------------------------+ *
#******************************************************************************
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.