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