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]]