Tanggal :September 25, 2020

R35. Fitting Data using Newton’s method in R

Spread the love

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 yn')
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:

Share

Leave a Reply

Your email address will not be published. Required fields are marked *