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]]
Hi r-users,
I have this code below but I don't understand the error message:
cumdensity <- 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
cden <- c1*t1bes1*t2
cden
}
pdensity <- function(z)
{ alp <- 2.0165; rho <- 0.868; a <- alp-0.5;
k1 <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)
k2 <- sqrt(rho)/(1-rho)
t1 <- exp(-z/(1-rho))
t2 <- (z/(2*k2))^a
bes1 <- besselI(z*k2,a)
bes2 <- besselI(z*k2,alp+0.5)
pp <- k1*t1*t2*(bes1/(rho-1) + a*bes1/z + k2*bes2 + a*bes1/z)
pp
}
## Newton iteration
newton_gam <- function(z)
{ n <- length(z)
r <- runif(n)
tol <- 1E-6
cdf <- vector(length=n, mode="numeric")
fprime <- vector(length=n, mode="numeric")
f <- vector(length=n, mode="numeric")
## Newton loop
for (i in 1:1000)
{ cdf <- cumdensity(z)
fprime <- pdensity(z)
f <- cdf - r
# Newton method
z <- z - f/fprime
if (f < tol) break
}
cbind(z,cdf)
}
alp <- 2.0165
bt1 <- 29.107 ; bt2 <- 41.517
r1 <- rgamma(10, shape=alp, scale = bt1)
r2 <- rgamma(10, shape=alp, scale = bt2)
z <- (r1/bt1)+(r2/bt2); sort(z)
OUTPUT
> z <- (r1/bt1)+(r2/bt2); sort(z)
[1] 0.9344932 1.3117707 1.4514991 2.2637735 2.2795669 3.0548586 3.4485707
4.1837583 4.2139719 4.5334232
> newton_gam(z)
z cdf
[1,] -25.540473 0.1883006
[2,] -2.079104 0.1725552
[3,] -2.878791 0.1331209
[4,] -81.317382 0.1881811
[5,] 8.395880 0.1376294
[6,] -20.381830 0.1601000
[7,] 7.124515 0.1682288
[8,] 16.099559 0.1755200
[9,] -2.635536 0.1218351
[10,] 6.438700 0.1341994
Warning message:
In if (f < tol) break :
the condition has length > 1 and only the first element will be used
What does this mean?
Thank you for any help given.
[[alternative HTML version deleted]]