[R] newton.method

[Berend Hasselman]

sammyny wrote:
> 
> I have this function:
> function(x)
> -0.3*x*exp(-(((log(x)+(0.03+0.3*0.3/2)*0.5)/(0.3*sqrt(0.5)))^2)/2)/(2*sqrt(2*pi*0.5))
> +
> 0.03*exp(-0.03*0.5)*pnorm(-(log(x)+(0.03-0.3*0.3/2)*0.5)/(0.3*sqrt(0.5)))
> ......
> 
> The newton method is not generating good result.
>> for( x in seq(0,8,1)) { z <- nleqslv(x,f); print(c(z$x,z$fvec))} 
> [1]           NaN 1.340781e+154
> [1]           NaN 1.340781e+154
> [1]  3.406704e+00 -5.608658e-09
> [1]  3.340717e+00 -9.480243e-09
> [1]  4.000000e+00 -5.470677e-11
> [1]  5.000000e+00 -3.386274e-14
> [1]  6.000000e+00 -3.559448e-17
> [1]  7.000000e+00 -6.065006e-20
> [1]  8.000000e+00 -1.581856e-22
> 
> In this case I am looking for largest value of 'x' such that f(x) >= 0
> 

The newton method (or rather Broyden to be more precise as you didn't use
method="Newton")
is finding values for x which make f(x) close to zero.

But you should follow the advice of Ravi and Dennis (several replies before
this reply) and plot the function.
I did that as follows:

curve(f, 0, 2)         
abline(0, 0, lty = 2)

and then all pieces of the puzzle fall into place.

-  Your function  goes to zero from x=2 onwards.
-  Your function has a zero or almost zero gradient for x in
c(0.0000001,0.5) 
-  Your function has a very steep gradient for x in c(0.5,1.0)

Whenever you starting value is slightly larger than 1.0 the initial
derivative of the function is > 0
>From the plot you can deduce that a Newton method will move away from the
solution you are looking for,
which is exactly what is happening for your starting values >= 2.

>From the plot you can also see that for a Newton method to succeed the
starting value should be between
0.6 (approximately) and .9.
For starting values x=0.6 and x=0.8 (look at the plot) nleqslv finds the
correct solution.

Summarizing, try this

for(x in seq(0.6,0.9,length=5)) { 
    cat(sprintf("Starting value for x=%g and f(x)=%g\n",x,f(x)))
    z <- nleqslv(x,f)
    print(c(z$x,z$fvec))                            
    cat("=========================================\n")
}

Newton/Broyden (or any method for that matter) need reasonably sensible
starting values otherwise they will fail. And you need to analyze/plot/.....
your function to find out the range for sensible starting values.

Berend
-- 
View this message in context: http://r.789695.n4.nabble.com/newton-method-tp2306111p2310422.html
Sent from the R help mailing list archive at Nabble.com.