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