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

