Nothing
snewtonm<-function(par,fn,gr,hess,control=list(trace=1, maxit=500),...) {
## Safeguarded Newton minimizer
##
##Input
## - fn is the function we wish to minimize
## - gr is the gradient function
## - hess is the hessian function
## - ... are exogenous data
##
##?? fixup documentation here??
## - par is the initial value
## - ... is data used in the function fn
##Output (list) -- need to match optim() output!! ???
## - xs is the value at the minimum
## - fv is the fn evaluated at xs
## - grd is the gradient
## - Hess is the Hessian
## - niter is the number of interations needed (gradient and Hessian evals).
## - add fevals??, other reports
npar <- length(par)
ctrldefault <- list(
trace = 0,
maxit = 500,
maxfevals = npar*500,
acctol = 0.0001,
epstol = .Machine$double.eps,
stepdec = 0.2,
stepmax = 5,
stepmin = 0,
offset = 100.0,
defstep=1,
bigval = .Machine$double.xmax*0.01,
watch = FALSE
)
ncontrol <- names(control)
nctrld <- names(ctrldefault)
for (onename in nctrld) {
if (! (onename %in% ncontrol)) {
control[onename]<-ctrldefault[onename]
}
}
trace <- control$trace # convenience
cat("trace =",trace,"\n")
cat(" Start snewtonm ")
nfn <- 0
ngr <- 0
nhess <- 0
eps0<-.Machine$double.eps
eps <- 10*eps0
fval <- fn(par, ...)
nfn <- nfn + 1
cat(" f0=",fval," at ")
print(par)
fbest <- abs(fval)*1.1 + 100.
itn <- 1
lambdamin<-(eps0^(1/4)) ## ?? do better
laminc <- 10.
lamdec <- 0.15
lambda<- lambdamin/lamdec## ?? do better
while ( (itn <= control$maxit) && (fval < fbest) ) { ## main loop
fbest <- fval
lambda <- lambda * lamdec
if (control$trace) {cat("Iteration ",itn," fbest=",fbest,"\n")}
grd<-gr(par,...)
ngr <- ngr + 1
if ( max(abs(grd)) < eps ) {
cat("Small gradient\n")
break
}
H<-hess(par,...)
nhess <- nhess + 1
while (fval >= fbest){
itn <- itn+1
Haug<-H + diag(npar)*lambda # To avoid singularity
stp<-solve(Haug, -grd)
xn <- par + stp # try unit step
if (identical(par,xn)) {
# if (max(abs(stp)) <= 1000*eps*max(abs(par))){
cat("Null step\n")
break
}
# if (control$trace) {cat(" step =", gvl," fval=", attr(gvl,"Fval"),"\n")}
fval <- fn(xn, ...)
nfn <- nfn + 1
if (control$trace) {cat(" lambda =", lambda," fval=", fval,"\n")}
if (fval >= fbest) {lambda <- max(lambdamin,lambda)*laminc} # increase lambda
}
# if (control$trace) {cat(" Success at lambda =", lambda," fval=", fval,"\n")}
par<-xn
if (control$watch) { tmp <- readline("end iteration") }
}
if (itn==control$maxit) {
print("NewtonR: Failed to converge!")
}
cat("Finished\n")
out<-NULL
out$xs<-xn
out$fv<-fn(xn,...)
out$grd<-grd
out$Hess<-H
out$counts <- list(niter=itn, nfn=nfn, ngr=ngr, nhess=nhess)
out
}
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.