# Estimate the West et al. model by gradient descent on nonlinear least squares
# Inputs: vector of initial parameter guesses, scale factor first, then
  # scaling exponent (params); vector of city populations (N); vector of city
  # per-capita economic outputs (Y); size of derivative of mean squared error
  # which is effectively zero (stopping.deriv); maximum number of iterations
  # (max.iterations); factor to multiply by the gradient (step.scale);
  # size of steps for taking numerical derivatives (deriv.increments)
# Outputs: List containing final value of the parameters, final value of
  # the gradient, the number of iterations used, and a flag for whether gradient
  # descent converged
my.nls.2 <- function(params, N=gmp$pop, Y=gmp$pcgmp, stopping.deriv,
  max.iterations, step.scale,deriv.increments) {
  mse <- function(params) {return(mean((Y-params[1]*N^params[2])^2)) }
  mse.grad <- function(params,deriv.increments) {
    p <- length(params); stopifnot(p==length(deriv.increments))
    gradient <- vector(length=p)
    mse.0 <- mse(params)
    for (i in 1:p) {
      new.params <- params
      new.params[i] <- params[i] + deriv.increments[i]
      new.mse <- mse(new.params)
      gradient[i] <- (new.mse-mse.0)/deriv.increments[i]
    }
    return(gradient)
  }
  for (iteration in 1:max.iterations) {
    gradient <- mse.grad(params,deriv.increments)
    params <- params - step.scale*gradient
    if(all(abs(gradient) < stopping.deriv)) { break() }
  }
  fit <- list(params=params,gradient=gradient,iterations=iteration,
    converged=(iteration < max.iterations))
  return(fit)
}