[R] Newton method again

[Berend Hasselman]

Roslina Zakaria wrote:
> 
> Hi Ravi,
> I did ask you some question regarding newton method sometime ago..  Now I
> have fixed the problem and I also wrote 2 looping code (ff1 and ff2) to
> evaluate the modified Bessel function of the first kind and call them in
> the newton code.  But I dont't understand why it gives the error message
> but still give the result but it is incorrect(too big and too small).
>  
> ff1 <- function(bb,eta,z){
> r <- length(z)
> for (i in 1:r) {
> sm <- sum(besselI(z[i]*bb,eta)/besselI(z[i]*bb,eta+1))*z[i]}
> sm
> }
> ff1(bb,eta,z)
> ff2 <- function(bb,eta,z,k){
> r <- length(z)
> for (i in 1:r) {
> sm1 <-
> sum((z[i]*bb/2)*(psigamma((0:k)+eta+1,deriv=0)/(factorial(0:k)*gamma((0:k)+eta+1))))
> sm2 <- sum((besselI(z[i]*bb,eta)*log(z[i]*bb/2) -
> sm1)/besselI(z[i]*bb,eta))}
> sm2
> }
> ff2(bb,eta,z,10)
>  
> newton.input3 <- function(pars)
> {  ##  parameters to be approximated , note: eta <- alpha3-0.5
>    eta   <- pars[1]
>    bt3   <- pars[2]
>    bt4   <- pars[3]
>    rho   <- pars[4]
>    b1    <- (pars[2]-pars[3])^2+4*pars[2]*pars[3]*pars[4]
>    b2    <- sqrt(b1)
>    bb    <- b2/(2*pars[2]*pars[3]*(1-pars[4]))
>    bf2   <-
> (pars[3]+2*pars[2]*pars[4]-pars[2])/(2*pars[2]^2*(pars[4]-1)*b2)
>    bf3   <-
> (pars[2]+2*pars[3]*pars[4]-pars[3])/(2*pars[3]^2*(pars[4]-1)*b2)
>    bf4   <-
> (2*pars[2]*pars[3]*pars[4]+pars[2]^2+pars[3]^2)/(2*pars[2]*pars[3]*(pars[4]-1)^2*b2)
>    zsm   <- sum(z)
>    psigm <- psigamma(pars[1]+0.5,deriv=0) 
>    pdz   <- log(prod(z))
>    erh   <- (1+2*pars[1])*(pars[4]-1)
>    brh1  <- 2*pars[2]*pars[3]*pars[4]+pars[2]^2+pars[3]^2 
>    brh2  <- 2*pars[2]*pars[3]*(pars[4]-1)^2 
>    k     <- 1000
>    ## function
>    fn1a  <- pdz -r*(2*psigm + log(b1))/2
>    fn2a  <- (pars[2]*r*erh+zsm)/(2*pars[2]^2*(1-pars[4]))
>    fn3a  <- (pars[3]*r*erh+zsm)/(2*pars[3]^2*(1-pars[4]))
>    fn4a  <-
> (pars[2]*pars[3]*r*erh+(pars[2]+pars[3])*zsm)/(-2*pars[2]*pars[3]*(pars[4]-1)^2)
>    
>    ## function that involve modified Bessel function of 1st kind 
>    fn1b  <- ff2(bb,eta,z,k)
>    fn2b  <- bf2*ff1(bb,eta,z)
>    fn3b  <- bf3*ff1(bb,eta,z)
>    fn4b  <-  bf4*ff1(bb,eta,z)
>    
>    ##  final function
>    fn1   <- fn1a + fn1b
>    fn2   <- fn2a + fn2b
>    fn3   <- fn3a + fn3b
>    fn4   <-  fn4a + fn4b
>    fval  <- c(fn1,fn2,fn3,fn4)
>    ## output
>    list(fval=fval)
> }
> library(BB)
> start <- c(0.7,0.8,0.6,0.4)
> dfsane(pars=start,fn=newton.input3)
> newton.input3(start)
> 
>> library(BB)
>> start <- c(0.7,0.8,0.6,0.4)
>> dfsane(pars=start,fn=newton.input3)
> Error in dfsane(pars = start, fn = newton.input3) : element 1 is empty;
>    the part of the args list of 'length' being evaluated was:
>    (par)
>> newton.input3(start)
> $fval
> [1]   103.0642   452.5835   823.6637 -1484.3209
> There were 50 or more warnings (use warnings() to see the first 50)
>>
>  
> Here is my data:
>> z
>  [1]  4.2 11.2  0.8 20.4 16.6  3.8  1.2  4.0 10.8 10.2  6.6 25.6 18.2  4.6
> 15.0  1.2 12.0 25.4  6.4  1.6  4.8 10.0  3.0
> [24]  7.0  1.8 15.0  8.6 11.2  5.4  1.8 23.2 10.8 25.4  6.0  6.0  5.0  1.4
> 11.0  8.4  7.4  6.4  2.6  8.6 15.8
>>
> 
> 

You could also try package nleqslv which implements Newton and Broyden
methods for solving systems of equations.

I have tried to run your problem but you are not providing all the
information required.
Moreover your example contains errors: for example where are the arguments
defined in the call of ff1  on the line "ff1(bb,eta,z)"  right after the
definition of ff1?

Where is the variable "r" used in the lines calculating fn1a, fn2a etc. in
function newton.input3?
Is it the same as in ff1 and ff2? length(z)?

When I insert r<-length(z) in newton.input3() I get the results shown in
your post for $fval.

The warnings are being given by factorial(0:k):  In factorial(0:k) : value
out of range in 'gammafn'

Why are you assigning pars[1], pars[2] etc to scalars and then afterwards
not or hardly using them?

You code is inefficient since you are calling ff1 in newton.input3  three
times with exactly the same input.

I have tried to run your code in nleqslv but it appears to run very slowly
so I can't help you any further at this point in time.

What is the purpose of the loop in function ff1

    for (i in 1:r) { 
    sm <- sum(besselI(z[i]*bb,eta)/besselI(z[i]*bb,eta+1))*z[i]} 

(on returning from the function sm will contain the value obtained for i=r)
?

Given the presentation of your problem, I cannot make head or tail of what
you are trying to do so I can't help you any further.

Berend Hasselman



-- 
View this message in context: http://www.nabble.com/newton-method-tp22653758p23875467.html
Sent from the R help mailing list archive at Nabble.com.