Nothing
#******************************************************************************
#* +------------------------------------------------------------------------+ *
#* | Function 'nlmest.NM', M-estimate of a nonlinear function. | *
#* | Using Nelder Mead, derivative free. | *
#* | 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. | *
#* | vm: variance matrix for generalized minimization if NULL not | *
#* | generalized. | *
#* | | *
#* +------------------------------------------------------------------------+ *
#******************************************************************************
nlmest.NM<-function(formula,data,start=getInitial(formula,data),robfunc,
control=nlr.control(tolerance=1e-4, minlanda=1 / 2 ^ 25, maxiter= 100 * length(start),robscale=T)
,vm=NULL,rm=eiginv(t(chol(vm))),delta=NULL,...)
{
tolerance <- control$tolerance
maxiter <- control$maxiter
minlanda <- control$minlanda
trace <- control$trace
datalist<-as.list(data)
p <- formula$p
n <- length(data[[1]])
loc.start <- start
if(trace) plot(data[[formula$independent]],data[[formula$dependent]])
Fault2 <- Fault() ###### no error, automatically will be created by the Fault object
if (is.null(loc.start$sigma)){
dt1 <- c(data,loc.start)
ht0 <- eval(formula,dt1)
dv <- as.numeric(ht0$predictor) - as.numeric(ht0$response)
if(! is.null(vm)) dv <- rm %*% dv
loc.start[["sigma"]] <- mad(dv) #mscale(dv)
}
th <- loc.start ## with sigma
sigma <- loc.start[["sigma"]]
b<- unlist(loc.start[names(loc.start)!="sigma"]) ## without sigma
if(is.null(delta)){
delta<-b * .1
delta[b==0] <- 0.1
}
theta <- rep(b,p) ## without sigma
dim(theta) <- c(p,p)
theta <- t(theta)
diag(theta) <- diag(theta)+delta
theta <- rbind(b-delta,theta) ## p+1 rows of parameters, exclude sigma
th <- theta ## include sigma
colnames(theta) <- names(start)[names(start)!="sigma"]
th <- cbind(th, rep(loc.start[["sigma"]],p+1))
colnames(th) <- names(loc.start)
y <- rep(0,p+1)
iterhist <- NULL
ht <- list(NULL)
for (i in 1:(p+1)){
dataloss <- th[i,]
if(is.null(vm)) ht[[i]] <- robloss(formula,data,start=dataloss ,robfunc,control=control,...) ## th: with sigma
else ht[[i]] <- robloss.gn(formula,data,dataloss ,robfunc,rm,control=control,...)
if(is.Faultwarn(ht[[i]])){
ht[[i]]@FN <- "nlmest.NL"
return(ht[[i]])
}
y[i] <- as.numeric(ht[[i]]$htheta)
iterhist <- rbind(iterhist,c(iteration = 0,objfnc = as.numeric(ht[[i]]$htheta),dataloss,converge = 0))
if(trace) lines(data[[formula$independent]],as.numeric(ht[[i]]$fmod$predictor))
}
eol <- F
iterate <- 0
psum <- apply(theta,2,sum) ### sum or mean ????? corect later
dt1 <- c(data,psum)
ft <- eval(formula,dt1)
dv <- as.numeric(ft$predictor) - as.numeric(ft$response)
if(! is.null(vm)) dv <- rm %*% dv
sigma <- mad(dv)
yold <- max(y)
#///*************************************************
while (iterate <= maxiter) #///********* Start Iteration ****************
{
o <- order(y) # don't need a full sort, only smallest and two largest
theta <- theta[o,] # so could be done more efficiently
th <- th[o,]
ht<-ht[o]
y <- y[o]
ilo <- 1
ihi <- p+1
inhi <- p
rtol <- 2*abs(y[ihi]-y[ilo])/(abs(y[ihi])+abs(y[ilo])+1e-8)
rtol2 <- sqrt(sum(th[1,]-th[2,])^2) / sqrt(sum(th[1,])^2)
iterhist <- rbind(iterhist,c(iteration = iterate,objfnc = y[ilo],th[ilo,],converge = min(rtol,rtol2) ))
if (rtol < tolerance | rtol2 < tolerance){
Fault2 <- Fault()
#result=list(parameters = as.list(th[ilo,]), objfnc=ht[[ilo]], history=iterhist,Fault=Fault2)
break
}
if (iterate >= maxiter){
Fault2 <- Fault(FN=1,FF = "optim.NM")
#result=list(parameters = as.list(th[ilo,]), objfnc=ht[[ilo]], history=iterhist,Fault=Fault2)
break
}
iterate <- iterate+2 # new point chosen by reflecting the worst current through the plane
# of the others
z <- smptry2(th,y,psum,data=data,ihi,-1,ht,vm=vm,rm=rm,sigma=sigma,control=control,robfunc=robfunc,formula=formula,...)
if(is.Fault(z)) return(z)
if (z[[1]] <= y[ilo]) { # new point is best--try going further
z <- smptry2(z[[4]],z[[2]],z[[3]],data=data,ihi,2,ht,vm=vm,rm=rm,sigma=sigma,control=control,robfunc=robfunc,formula=formula,...)
yold<-y[1]
y <- z[[2]]; psum <- z[[3]]; th <- z[[4]];ht<-z[[5]]
sigma <- z$sigma
}
else
if (z[[1]] >= y[inhi]) {
ysave <- z[[2]][ihi] #new point is still worst, try smaller step
z <- smptry2(z[[4]],z[[2]],z[[3]],data=data,ihi,0.5,ht,vm=vm,rm=rm,sigma=sigma,control=control,robfunc=robfunc,formula=formula,...)
yold<-y[1]
y <- z[[2]]; psum <- z[[3]]; th <- z[[4]];ht<-z[[5]] # a final psum is that after it not any smptry compute
if (z[[1]] >= ysave) { # still bad, shrink simplex
for (i in (1:(p+1))[-ilo]) {
psum <- (theta[i,]+theta[ilo,])/2
theta[i,] <- psum
th[i,] <- c(psum,sigma)
if(is.null(vm)) ht[[i]] <- robloss(formula,data,start=th[i,] ,robfunc=robfunc,control=control,...) ## th: with sigma
else ht[[i]] <- robloss.gn(formula,data,th[i,] ,robfunc,rm,control=control,...)
if(is.Faultwarn(ht[[i]])){
ht[[i]]$Fault$FN <- "optim.NL"
return(ht[[i]])
}
y[i] <- as.numeric(ht[[i]]$htheta)
}
iterate <- iterate+p
psum <- apply(th[,1:p],2,sum)
dt1 <- c(data,psum)
ft <- eval(formula,dt1)
dv <- as.numeric(ft$predictor) - as.numeric(ft$response)
if(! is.null(vm)) dv <- rm %*% dv
sigma <- mscale(dv)
}
else sigma <- z$sigma # anytime the final psum must compute sigma
}
else {
y <- z[[2]]; psum <- z[[3]]; th <- z[[4]];ht<-z[[5]]
sigma <- z$sigma
iterate <- iterate-1
}
if(trace) lines(data[[formula$independent]],ht[[ilo]]$fmod$predictor)
} #### end itteration #####
##################################
#cat("finished iteration=",iterate ,"\n")
#//////////*******************************************************************************************************************
#//////////*******************************************************************************************************************
#//////////******** preparing output ***********************************************************************************
#//////////******** ***********************************************************************************
#### result=list(parameters = as.list(th[ilo,]), objfnc=ht[[ilo]], history=iterhist,Fault=Fault2)
yresp <- as.numeric(ht[[ilo]]$fmod$response)
ypred <- as.numeric(ht[[ilo]]$fmod$predictor)
ybar <- mean(yresp)
nlrho <- 1 - sum( ( yresp - ypred ) ^ 2 ) /
sum( (ypred - ybar ) ^ 2 )
# curv1 <- curvature(gradient = attr(ht[[ilo]]$fmod$predictor,"gradient"),
# hessian = attr(ht$fmod$predictor,"hessian"),
# sigma = th[ilo,"sigma"])
if(is.null(vm))
result <- nl.fitt.rob(
parameters = as.list(th[ilo,]),
scale = sigma,
correlation = nlrho,
form = formula,
response = ht[[ilo]]$fmod$response,
predictor = ht[[ilo]]$fmod$predictor,
# curvature = curv1,
history = iterhist,
method = fittmethod(methodID= 7,
methodBR= 13, ### rob or nopnrob variance
detailBR= "Nelder Mead",
subroutine= "nlmest.NM"),
data = as.list(data),
sourcefnc = match.call(),
Fault = Fault2,
htheta = ht[[ilo]]$htheta,
rho = ht[[ilo]]$rho,
ri = ht[[ilo]]$ri,
curvrob = NULL,
robform = robfunc
)
else
result <- nl.fitt.rgn(
parameters = as.list(th[ilo,]),
scale = sigma,
correlation = nlrho,
form = formula,
response = ht[[ilo]]$fmod$response,
predictor = ht[[ilo]]$fmod$predictor,
# curvature = curv1,
history = iterhist,
method = fittmethod(methodID= 8,
methodBR= 13, ### rob or nopnrob variance
detailBR= "Nelder Mead",
subroutine= "nlmest.NM"),
data = as.list(data),
sourcefnc = match.call(),
Fault = Fault2,
htheta = ht[[ilo]]$htheta,
rho = ht[[ilo]]$rho,
ri = ht[[ilo]]$ri,
curvrob = NULL,
robform = robfunc,
vm = vm,
rm= rm,
gresponse = transformNR(ht[[ilo]]$fmod$response,rm),
gpredictor = transformNR(ht[[ilo]]$fmod$predictor,rm)
)
return(result)
}
#******************************************************************************
#* +------------------------------------------------------------------------+ *
#* | End of nlmest.NM' | *
#* | | *
#* | Hossein Riazoshams, Dep Stat, Stock | *
#* | | *
#* | Recoded from 'optim.NM' in March 20103 | *
#* | | *
#* +------------------------------------------------------------------------+ *
#******************************************************************************
##
##
###### function smptry ######
smptry2 <- function(th,y,psum,data,ihi,fac,ht,vm,rm,sigma,formula,control,...) { #### sigma is the scale computed at psum
ndim <- ncol(th)-1
n <- length(data[[1]])
fac1 <- (1-fac)/ndim
fac2 <- fac1-fac
ptry <- psum*fac1-th[ihi,1:ndim]*fac2
th2 <- as.list(ptry)
th2$sigma <- sigma
if(is.null(vm)) htobj <- robloss(formula,data,start=th2 ,...) ## th: with sigma
else htobj <- robloss.gn(formula,data,th2,rmat=rm,control=control,...)
if (is.Fault(htobj)) return(htobj)
tmp <- as.numeric(htobj$ri)
if(control$robscale) sigma<-mscale(tmp) # <- nl.mscale(tmp,robfunc,...)
else sigma <- sum(abs(tmp)) / n
ytry <- as.numeric( htobj$htheta)
if (ytry < y[ihi]) {
y[ihi] <- ytry
psum <- psum-th[ihi,1:ndim]+ptry
th[ihi,1:ndim] <- ptry
ht[[ihi]]<-htobj
dt1 <- c(data,psum)
ht0 <- eval(formula,dt1)
dv <- as.numeric(ht0$predictor) - as.numeric(ht0$response)
if(! is.null(vm)) dv <- rm %*% dv
sigma <- mscale(dv)
}
list(ytry,y,psum=psum,th=th,ht=ht,sigma=sigma)
}
###############################################
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.