R35. Fitting Data using Newton's method in R

Test vectors, x and y are created. x, the input, goes from 0 to 100 in steps of 1.

y, the output, corresponds to the function x/(25+x) with some added noise.

Now when using Newton's method, we have a vector of length 101. Some of the values will be positive and some negative. The function fit finds the sum over all the terms. Finding the root of fit, will find the curve in the middle of the data points. We find the best fit according to equation x/(k+x). We know that k has to be near 25. There are many functions for finding best fits in R.

In many cases, Newton's method will not converge such as when the y' keeps on changing sign. In those cases you can always plot the function over a grid of points and then find an approximate local minima, and use that as the starting value in Newton's method (the argument x0).

# ex35.R
fit <- function(k) {
  sum(y-x/(k+x))
}

newton <- function(f, tol=1E-10, x0=1, N=50, h=.001) {
  for (i in 1:N) {
    df.dx <- (f(x0+h/2)-f(x0-h/2))/h
    x1 <- x0 - f(x0)/df.dx
    error <- abs(x1-x0)
    cat("i = ",i,", x1 = ",x1,", error = ",error,'\n')
    if (error < tol) break
    x0 <- x1
  }
  x1
}

cat('Test data x and y\n')
x <- 0:100
y <- x/(25+x)+ runif(length(x), -.1, .1)
plot(x, y, type='l',
     col = 'blue',
     main = expression("Fitting y = x/(k"[0]*"+x)"))
k0 <- newton(fit)
cat('k0 =',k0,'\n')
legend("center",  pch=c('--'),
        col = c("blue", "magenta"), 
        legend = c("actual","fit"))
lines(x, x/(k0+x),col='magenta')

#Test data x and y
#i =  1 , x1 =  11.33843 , error =  10.33843 
#i =  2 , x1 =  22.53979 , error =  11.20136 
#i =  3 , x1 =  26.25 , error =  3.710206 
#i =  4 , x1 =  26.49259 , error =  0.2425991 
#i =  5 , x1 =  26.4935 , error =  0.0009042225 
#i =  6 , x1 =  26.4935 , error =  1.245337e-08 
#i =  7 , x1 =  26.4935 , error =  1.065814e-14 
#k0 = 26.4935 

Output: