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