Indexing fix for small sigma

R/BB.R

@@ -85,7 +85,11 @@ dBB <- function(x, mu = 0.5, sigma = 1, bd = 10, log = FALSE)
                   +lgamma((1/sigma))+lgamma(xx+mu*(1/sigma))
                   +lgamma(bd+((1-mu)/sigma)-xx)-lgamma(mu*(1/sigma))
                   -lgamma((1-mu)/sigma)-lgamma(bd+(1/sigma)))
-logfy[sigma<0.0001]  <- dBI(xx, mu = mu, bd=bd, log = TRUE)
+       idx <- sigma < 0.0001
+       if (any(idx)){
+        logfy[idx]  <- dBI(xx[idx], mu = mu[idx], bd=bd[idx], log = TRUE)
+       }
+
          fy <- if(log == FALSE) exp(logfy) else logfy
   fy[x < 0] <- 0 
 fy[x > bd] <- 0 

R/DELAPORT.R

@@ -160,7 +160,11 @@ logpy0 <- -mu*nu-(1/sigma)*(log(1+mu*sigma*(1-nu)))
  logfy <-  logpy0-lgamma(x+1)+S
  if(log==FALSE) fy <- exp(logfy) else fy <- logfy
    fy[sigma>0.0001] <- fy
-   fy[sigma<=0.0001] <- dPO(x, mu = mu, log = log) 
+   idx <- sigma <= 0.0001
+   if (any(idx)){
+    fy[idx] <- dPO(x[idx], mu = mu[idx], log = log)
+   }
+
    fy[x < 0] <- 0
    fy[x == Inf] <- 0
   fy

R/GeneralisedPoisson.R

@@ -90,7 +90,11 @@ dGPO<-function(x, mu = 1, sigma = 1, log = FALSE)
      sigma <- rep(sigma, length = ly)
         mu <- rep(mu, length = ly)   
       logL <- x*log(mu/(1+sigma*mu))+(x-1)*log(1+sigma*x)+(-mu*(1+sigma*x))/(1+sigma*mu)-lgamma(x+1)
-logL[sigma>0.000001]  <- dPO(x, mu = mu, log = TRUE)   # if sigma too small Poisson 
+      idx <- sigma < 0.000001
+      if (any(idx)){
+       logL[idx]  <- dPO(x[idx], mu = mu[idx], log = TRUE)   # if sigma too small Poisson 
+      }
+
           Lik <- if (log) logL else exp(logL)
    Lik[x < 0] <- 0
 Lik[x >= Inf] <- 0
@@ -123,7 +127,11 @@ y.y[qq[i]==Inf] <- 1
          FFF[i] <- sum(pdfall)                                             
     }  
     cdf <- FFF
-cdf[sigma<0.0001]  <- pPO(q, mu = mu, log.p = TRUE)
+    idx <- sigma < 0.0001
+    if (any(idx)){
+     cdf[idx]  <- pPO(q[idx], mu = mu[idx], log.p = TRUE)   
+    }
+
        cdf[q<0] <- 0
     cdf[q>=Inf] <- 1 
             cdf <- if(lower.tail==TRUE) cdf else 1-cdf

R/NBF.R

@@ -119,7 +119,11 @@ dNBF<-function(x, mu=1, sigma=1, nu=2, log=FALSE)
      mu1 <- mu
   sigma1 <- sigma*mu^(nu-2)
       fy <-  dnbinom(x, size=1/sigma1, mu = mu1, log = log)
-fy[sigma1<=0.0001] <-  dPO(x, mu = mu1, log = log) 
+      idx <- sigma1 <= 0.0001
+      if (any(idx)){
+       fy[idx] <-  dPO(x[idx], mu = mu1[idx], log = log) 
+      }
+
   fy[x < 0] <- 0 
   fy[x == Inf] <- 0 
   fy
@@ -140,8 +144,15 @@ pNBF <- function(q, mu=1, sigma=1, nu=2, lower.tail = TRUE, log.p = FALSE)
      cdf <- rep(0,length = ly) 
      mu1 <- mu
   sigma1 <- sigma*mu^(nu-2)
- cdf[sigma1>0.0001] <- pnbinom(q, size=1/sigma1, mu=mu1, lower.tail=lower.tail,log.p=log.p)
-cdf[sigma1<=0.0001] <- ppois(q, lambda = mu1, lower.tail = lower.tail, log.p = log.p) 
+  idx <- sigma1 > 0.0001
+  if (any(idx)){
+   cdf[idx] <- pnbinom(q[idx], size=1/sigma1[idx], mu=mu1[idx], lower.tail=lower.tail,log.p=log.p) 
+  } 
+  idx <- sigma1 <= 0.0001
+  if (any(idx)){
+    cdf[idx] <- ppois(q[idx], lambda = mu1[idx], lower.tail = lower.tail, log.p = log.p)
+  }
+
   cdf[q < 0] <- 0 
   cdf[q >= Inf] <- 1
   cdf

R/NBI.R

@@ -68,8 +68,12 @@ if (any(sigma <= 0) )  stop(paste("sigma must be greater than 0 ", "\n", ""))
                mu <- rep_len(mu, n)
             sigma <- rep_len(sigma, n)
                fy <- rep_len(0,n)        
-               fy <-  dnbinom(x, size=1/sigma, mu = mu, log = log) 
-fy[sigma<=0.0001] <-  dPO(x, mu = mu, log = log)
+               fy <-  dnbinom(x, size=1/sigma, mu = mu, log = log)
+               idx <- sigma <= 0.0001
+               if (any(idx)){
+                fy[idx] <-  dPO(x[idx], mu = mu[idx], log = log)
+               } 
+
        fy[x < 0]  <- 0 
     fy[x == Inf]  <- 0 
         return(fy)
@@ -90,8 +94,12 @@ if (any(sigma <= 0) )  stop(paste("sigma must be greater than 0 ", "\n", ""))
               cdf <- rep_len(0,n)
              cdf  <-  pnbinom(q, size=1/sigma, mu=mu, 
                               lower.tail=lower.tail, log.p = log.p)
-cdf[sigma<=0.0001]<-  ppois(q, lambda = mu, 
-                            lower.tail = lower.tail, log.p = log.p)            
+             idx <- sigma <= 0.0001
+             if (any(idx)){
+              cdf[idx]<-  ppois(q[idx], lambda = mu[idx], 
+                            lower.tail = lower.tail, log.p = log.p)  
+             }
+
         cdf[q < 0] <- 0
         cdf[q == Inf] <- 1
         return(cdf)

R/NBII.R

@@ -93,7 +93,11 @@ if (any(sigma <= 0) )  stop(paste("sigma must be greater than 0 ", "\n", ""))
             sigma <- rep_len(sigma, n)
                fy <- rep_len(0,n)        
                fy <- dnbinom(x, size=mu/sigma, mu=mu, log=log) 
-fy[sigma<=0.0001] <-  dpois(x, lambda = mu, log = log)
+               idx <- sigma <= 0.0001
+               if (any(idx)){
+                fy[idx] <-  dpois(x[idx], lambda = mu[idx], log = log)
+               }
+
         fy[x < 0] <- 0 
      fy[x == Inf] <- 0 
   return(fy)