Add probability mass functions

DESCRIPTION

@@ -18,7 +18,9 @@ LinkingTo:
     Rcpp
 Imports: 
     Rcpp,
-    mcglm (>= 0.5.0)
+    mcglm (>= 0.5.0),
+    statmod,
+    tweedie
 RoxygenNote: 6.1.1
-Remotes:  
+Remotes: 
     wbonat/mcglm

NAMESPACE

@@ -2,8 +2,15 @@
 
 export(compute_logk)
 export(compute_logz)
+export(dcmp)
+export(ddpo)
+export(ddwe)
+export(dgct)
+export(dgpo)
+export(dptw)
 export(flexcm)
 importFrom(Rcpp,sourceCpp)
+importFrom(stats,dpois)
 importFrom(stats,glm.fit)
 importFrom(stats,model.frame)
 importFrom(stats,model.matrix)

R/probability.R

@@ -0,0 +1,134 @@
+#-----------------------------------------------------------------------
+#' @title COM-Poisson probability mass function
+#' @param y a positive integer value.
+#' @param mu the (approximate) mean parameter.
+#' @param nu the dispersion parameter.
+#' @param log logical. If \code{TRUE} it returns the logarithm of pmf.
+#' @author Eduardo Jr <[email protected]>
+#' @export
+#'
+dcmp <- function(y, mu, nu, log = FALSE) {
+  loglambda <- suppressWarnings(
+    nu * log(mu + (nu - 1) / (2 * nu))
+  )
+  # Get the normalizing constants
+  logz <- compute_logz(loglambda, mu, nu)
+  # Compute the loglikelihood
+  pmf <- y * loglambda - nu * lfactorial(y) - logz
+  if (!log) pmf <- exp(pmf)
+  return(pmf)
+}
+
+#-----------------------------------------------------------------------
+#' @title Gamma-count probability mass function
+#' @param y a positive integer value.
+#' @param kappa the (asymptote) mean parameter.
+#' @param alpha the dispersion parameter.
+#' @param log logical. If \code{TRUE} it returns the logarithm of pmf.
+#' @author Eduardo Jr <[email protected]>
+#' @export
+#'
+dgct <- function (y, kappa, alpha, log = FALSE) {
+  alpha_kappa <- alpha * kappa
+  alpha_y <- alpha * y
+  pmf <- pgamma(1L, shape = alpha_y, rate = alpha_kappa) -
+    pgamma(1L, shape = alpha_y + alpha, rate = alpha_kappa)
+  if (log) pmf <- log(pmf)
+  return(pmf)
+}
+
+#-----------------------------------------------------------------------
+#' @title Discrete Weibull probability mass function
+#' @param y a positive integer value.
+#' @param q the "location" parameter (\code{0<q<1}).
+#' @param rho the dispersion parameter.
+#' @param log logical. If \code{TRUE} it returns the logarithm of pmf.
+#' @author Eduardo Jr <[email protected]>
+#' @export
+#'
+ddwe <- function(y, q, rho, log = FALSE) {
+  logq <- log(q)
+  pmf <- exp(y^rho * logq) - exp((y + 1)^rho * logq)
+  if (log) pmf <- log(pmf)
+  return(pmf)
+}
+
+#-----------------------------------------------------------------------
+#' @title Generalized Poisson probability mass function
+#' @param y a positive integer value.
+#' @param mu the mean parameter.
+#' @param sigma the dispersion parameter.
+#' @param log logical. If \code{TRUE} it returns the logarithm of pmf.
+#' @author Eduardo Jr <[email protected]>
+#' @export
+#'
+dgpo <- function (y, mu, sigma, log = FALSE) {
+  sigma_mu <- 1 + sigma * mu
+  sigma_y <- 1 + sigma * y
+  lsigma_mu <- suppressWarnings(log(sigma_mu))
+  lsigma_y <- suppressWarnings(log(sigma_y))
+  sy_sm <- suppressWarnings(sigma_y / sigma_mu)
+  pmf <- y * (log(mu) - lsigma_mu) +
+    (y - 1) * lsigma_y - mu *
+    (sy_sm) - lfactorial(y)
+  pmf[is.nan(pmf)] <- -Inf
+  if (!log) pmf <- exp(pmf)
+  return(pmf)
+}
+
+#-----------------------------------------------------------------------
+#' @title Double Poisson probability mass function
+#' @param y a positive integer value.
+#' @param mu the (approximate) mean parameter.
+#' @param phi the dispersion parameter.
+#' @param log logical. If \code{TRUE} it returns the logarithm of pmf.
+#' @author Eduardo Jr <[email protected]>
+#' @export
+#'
+ddpo <- function (y, mu, phi, log = FALSE) {
+  # Get the normalizing constants
+  logk <- compute_logk(mu, log(mu), phi, log(phi))
+  ly <- log(y); ly[y == 0] <- 1
+  pmf <- -0.5 * log(phi) - mu/phi - y + y * ly - lfactorial(y) +
+    (y/phi) * (1 + log(mu) - ly) - logk
+  if (!log) pmf <- exp(pmf)
+  return(pmf)
+}
+
+#-----------------------------------------------------------------------
+#' @title Poisson-Tweedie probability mass function
+#' @param y a positive integer value.
+#' @param mu the mean parameter.
+#' @param omega the dispersion parameter.
+#' @param power the power parameter.
+#' @param npts number of points in the gaussian quadrature.
+#' @param log logical. If \code{TRUE} it returns the logarithm of pmf.
+#' @author Eduardo Jr <[email protected]>
+#' @importFrom stats dpois
+#' @export
+#'
+dptw <- function(y, mu, omega, power, npts = 100L, log = FALSE) {
+    pts <- statmod::gauss.quad(npts, kind = "laguerre")
+    pjoint <- function(y, z) {
+        dpois(y, lambda = z) *
+            tweedie::dtweedie(z, mu = mu, phi = omega, power = power)
+    }
+    pmarginal <- function(y) {
+        vapply(y, FUN = function(yi) {
+            fpts <- pjoint(y = yi, z = pts$nodes)
+            sum(pts$weights * fpts / exp(-pts$nodes))
+        }, FUN.VALUE = numeric(1))
+    }
+    out <- vapply(y, FUN = function(yi) {
+        if (yi == 0) {
+            se <- sqrt(mu + omega * mu^power)
+            yrange <- 1:min(1000, round(mu + 10 * se))
+            probs <- pmarginal(yrange)
+            out <- 1 - sum(probs)
+        } else {
+            pmarginal(yi)
+        }
+    }, FUN.VALUE = numeric(1))
+    if (log) out <- log(out)
+    return(out)
+}