Hi R-users,
 
I have this code here:

library(numDeriv)
 
fprime <- function(z)
{ alp  <- 2.0165;
  rho  <- 0.868;
 
# simplified expressions
  a      <- alp-0.5
  c1     <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)
  c2     <- sqrt(rho)/(1-rho)
  t1     <- exp(-z/(1-rho))
  t2     <- (z/(2*c2))^a
  bes1   <- besselI(z*c2,a)
  t1bes1 <- t1*bes1
  c1*t1bes1*t2
}
 
## Newton iteration
newton_gam <- function(z)
{ n   <- length(z)
  r   <- runif(n)
  tol <- 1E-6
  cdf <- vector(length=n, mode="numeric")
 
  for (i in 1:1000)
  { # numerical intergration to find the cdf
    cdf  <- integrate(fprime, lower = 0, upper = z)$value
 
    # Newton method
    z    <- z - (cdf - r)/fprime(z)
 
    if (tol < 1e-10) break
   }
  cbind(z,cdf)
}
 
bt1 <- 29.107 ; bt2 <- 41.517
x1  <- 30; x2 <- 10;
 
z <- (x1/bt1)+(x2/bt2); z
newton_gam(z)
 
> bt1 <- 29.107 ; bt2 <- 41.517
> x1  <- 30; x2 <- 10;
> z <- (x1/bt1)+(x2/bt2); z
[1] 1.271545
> newton_gam(z)
            z       cdf
[1,] 4.128138 0.6065616
 
But when I try with different values of X1 and X2, it gives
 
> x1  <- 1; x2 <- 5;
> z <- (x1/bt1)+(x2/bt2); z
[1] 0.1547886
> newton_gam(z)
Error in integrate(fprime, lower = 0, upper = z) : 
  non-finite function value
In addition: There were 22 warnings (use warnings() to see them)
 
warnings()
Warning messages:
1: In besselI(z * c2, a) : value out of range in 'bessel_i'
2: In besselI(z * c2, a) : value out of range in 'bessel_i'
3: In besselI(z * c2, a) : value out of range in 'bessel_i'
4: In besselI(z * c2, a) : value out of range in 'bessel_i'
5: In besselI(z * c2, a) : value out of range in 'bessel_i'
 
I try checking the besselI(z * c2, a) values, it seem that no problem.  Why it gives such warning? Any idea?
 
> x1  <- 0.5; x2 <- 0.8
 
> c2 <- sqrt(rho)/(1-rho)
 
> z  <- (x1/bt1)+(x2/bt2)
 
> besselI(z*c2,a)
[1] 0.03337589
> source(.trPaths[5], echo=TRUE, max.deparse.length=10000)
 
> x1  <- 1; x2 <- 5
 
> c2 <- sqrt(rho)/(1-rho)
 
> z  <- (x1/bt1)+(x2/bt2)
 
> besselI(z*c2,a)
[1] 0.3339742
 
> x1  <- 10; x2 <- 50
 
> c2 <- sqrt(rho)/(1-rho)
 
> z  <- (x1/bt1)+(x2/bt2)
 
> besselI(z*c2,a)
[1] 6076.817
 
> x1  <- 100; x2 <- 200
 
> c2 <- sqrt(rho)/(1-rho)
 
> z  <- (x1/bt1)+(x2/bt2)
 
> besselI(z*c2,a)
[1] 1.018641e+24
 
Thank you.
 


      
	[[alternative HTML version deleted]]

