hwk 8.4

gistfile1.txt

data <- read.table("~/Dropbox/numericalmethods/nonlinear.txt", header=TRUE, quote="\"")

#param order: d, r
fxn.4 <- list(
  f.x = function(params, data) return( sum((data$y - params[1]*exp(-params[2] * data$x) )^2) ),
  grad = function(params, data, f.x, h) { 
    d.d <- (f.x(params = c(params[1] + h, params[2]), data) - f.x(params, data)) / h
    d.r <- (f.x(params = c(params[1], params[2] + h), data) - f.x(params, data)) / h
    return(c(d.d,d.r))
    },
  hess = function(params, data, f.x) {
    return( hessian(f.x, x=params, data=data) ) 
    },
  norm = function(x) return(sqrt(sum(x)^2)) ,
  graph = function(min, data) {
  qplot(x = data$x, y = data$y) + 
  stat_function(fun = function(x) min[1]*exp(-min[2] * x), geom="line", lwd = 1, aes(colour="(0,0)")) +
  xlab("X") + ylab("Y")
  }
    )

newton <- function(start, err, fxn=fxn.4, h = 10^-3, back = TRUE, grapher = FALSE, data=data, adj = .001) {
  out <- list("min" = NULL, "steps" = 0, "stepsize" = NULL, "start" = start, "path" = start, "graph" = NULL)
  x.i <- out$start
  
  while(fxn$norm(fxn$grad(x.i, data, fxn$f.x, h)) > err) {
    step <- 1
    M <- fxn$hess(x.i, data, fxn$f.x)
    
    #checking to see if any eigenvalues are negative, and if negative, shift by some lambda. 
    lambda <- min(eigen(M)$values)
    if(lambda < 0) M <- (abs(lambda) + adj)*diag(rep(1,nrow(M))) + M 
    
    #calculate direction with multiple dimensions
    dir <- -solve(M)%*%fxn$grad(x.i, data, fxn$f.x, h)
    
    #backtracking, can be turned on and off
    if(back) while(fxn$f.x((x.i + step*dir), data) > fxn$f.x(x.i, data)) step = step/2
    
    x.i <- x.i + step*dir
    out$path <- rbind(out$path, as.numeric(x.i) )
    out$steps <- out$steps + 1
    out$stepsize <- rbind(out$stepsize, step)
  }
  out$min <- as.numeric(x.i)

  if(grapher) out$graph <- fxn$graph(out$min, data) #graph path of results
  
  return(out)
}

newton(start = c(10,10), err = .001, fxn=fxn.4, h = 10^-4, back = TRUE, grapher = FALSE, data=data, adj = .001)