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
>
 
The answer that I should get is c(0.4960, 6.6265,2.6470,0.4378)
 
Thank you so much for help and time.


      
	[[alternative HTML version deleted]]


