Hi r-users,
 
I hope somebody can help me with this code. 
I would like to solve for z values using newton iteration method.  I 'm not 
sure which part of the code is wrong since I'm not very good at programming but 
would like to learn.  There seem to be some output but what I expected is a 
vector of z values.  Thank you so much for any help given.
 
newton.inputsingle <- function(pars,n)
{  runi    <- runif(974, min=0, max=1)
   lendt   <- length(runi)
   ## Parameter to estimate
   z <- vector(length=lendt, mode= "numeric")
   z  <- pars[1]
   
   ## Constant value  
   
   alp  <- 2.0165 ; rho <- 0.868; 
   c    <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)
   
   for (i in 1:n)
   {  t1   <- exp(-pars[1]/(1-rho))                       
      t2   <- (pars[1]*(1-rho)/(2*sqrt(rho)))^(alp-0.5)   
      bes1 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp-0.5)  
      bes2 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp-1.5)
      bes3 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp+0.5)
 
   ## Equation
      f   <- c*t1*t2*bes1 - runi
   
   ## derivative
      fprime   <- c*t1*t2*( -bes1/(1-rho) + (alp-0.5)*bes1/pars[1] + 
sqrt(rho)*(bes2-bes3)/(2*(1-rho)))
      z[i+1] <- z[i] - f/fprime 
      }
      z
}
 
pars <- 0.5          
newton.inputsingle(pars,5)
 
The output :
 
> pars <- 0.5          
> newton.inputsingle(pars,5)
[1]  0.5000000 -0.4826946 -1.4653892 -2.4480838 -3.4307784 -4.4134730
Warning messages:
1: In z[i + 1] <- z[i] - f/fprime :
  number of items to replace is not a multiple of replacement length
2: In z[i + 1] <- z[i] - f/fprime :
  number of items to replace is not a multiple of replacement length
3: In z[i + 1] <- z[i] - f/fprime :
  number of items to replace is not a multiple of replacement length
4: In z[i + 1] <- z[i] - f/fprime :
  number of items to replace is not a multiple of replacement length
5: In z[i + 1] <- z[i] - f/fprime :
  number of items to replace is not a multiple of replacement length




      
        [[alternative HTML version deleted]]

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to