Init the package (with C++ functions)

.Rbuildignore

@@ -0,0 +1,3 @@
+^data-raw$
+^LICENSE\.md$
+^README\.md$

.gitignore

@@ -0,0 +1,15 @@
+# Backup/History files
+*~
+.#*
+#*
+.Rhistory
+
+# Rstudio projects
+.Rproj.user
+*.Rproj
+
+# Others
+*.Rout
+*.html
+roteiro*.R
+inst/doc

DESCRIPTION

@@ -0,0 +1,21 @@
+Package: flexcm
+Title: Flexible Models for Dispersed Count Data
+Version: 0.0.1.9000
+Authors@R:
+    person(given = "Eduardo",
+           family = "Elias Ribeiro Junior",
+           role = c("aut", "cre"),
+           email = "[email protected]")
+Description: Provide functions to fit flexible count models that can
+  handle equi-, over-, and underdispersion, namely we considered
+  COM-Poisson, Gamma-count, discrete Weibull, generalized Poisson,
+  double Poisson and Poisson-Tweedie models. The normalizing constants
+  for COM-Poisson and double Poisson models are written in C++.
+License: GPL-3
+Encoding: UTF-8
+LazyData: true
+LinkingTo: 
+    Rcpp
+Imports: 
+    Rcpp
+RoxygenNote: 6.1.1

NAMESPACE

@@ -0,0 +1,5 @@
+# Generated by roxygen2: do not edit by hand
+
+export(compute_logk)
+export(compute_logz)
+useDynLib(flexcm, .registration = TRUE)

R/RcppExports.R

@@ -0,0 +1,36 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#' @title Compute COM-Poisson normalizing constante in log scale
+#' @description Computes the normalizing constant in the log scale
+#'   usign the LogSumExp trick to avoid numerical issues (See details).
+#' @details \code{logspace_add(a, b) = log(exp(a) + exp(b)) = b + log(1 +
+#'   exp(a - b))}, where \code{b > a}.
+#' @param loglambda A vector of logarithm of the \eqn{\lambda}
+#'   parameter.
+#' @param nu A vector of dispersion parameters \eqn{\nu}.
+#' @return The normalizing constant.
+#' @references Wikipedia, LogSumExp. \url{http://rstudio.com}.
+#' @author Eduardo Jr <[email protected]>
+#' @export
+compute_logz <- function(loglambda, nu) {
+    .Call(`_flexcm_compute_logz`, loglambda, nu)
+}
+
+#' @title Compute Double Poisson normalizing constante in log scale
+#' @description Computes the normalizing constant in the log scale
+#'   usign the LogSumExp trick to avoid numerical issues (See details).
+#' @details \code{logspace_add(a, b) = log(exp(a) + exp(b)) = b + log(1
+#'   + exp(a - b))}, where \code{b > a}.
+#' @param mu a vector of the \eqn{\mu} parameter.
+#' @param lmu logarithm of \eqn{\mu}.
+#' @param phi value of \eqn{\varphi} parameter.
+#' @param lphi logarithm of \eqn{\varphi}.
+#' @return The normalizing constant.
+#' @references Wikipedia, LogSumExp. \url{http://rstudio.com}.
+#' @author Eduardo Jr <[email protected]>
+#' @export
+compute_logk <- function(mu, lmu, phi, lphi) {
+    .Call(`_flexcm_compute_logk`, mu, lmu, phi, lphi)
+}
+

src/.gitignore

@@ -0,0 +1,3 @@
+*.o
+*.so
+*.dll

src/RcppExports.cpp

@@ -0,0 +1,44 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+// compute_logz
+NumericVector compute_logz(NumericVector loglambda, double nu);
+RcppExport SEXP _flexcm_compute_logz(SEXP loglambdaSEXP, SEXP nuSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type loglambda(loglambdaSEXP);
+    Rcpp::traits::input_parameter< double >::type nu(nuSEXP);
+    rcpp_result_gen = Rcpp::wrap(compute_logz(loglambda, nu));
+    return rcpp_result_gen;
+END_RCPP
+}
+// compute_logk
+NumericVector compute_logk(NumericVector mu, NumericVector lmu, double phi, double lphi);
+RcppExport SEXP _flexcm_compute_logk(SEXP muSEXP, SEXP lmuSEXP, SEXP phiSEXP, SEXP lphiSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type lmu(lmuSEXP);
+    Rcpp::traits::input_parameter< double >::type phi(phiSEXP);
+    Rcpp::traits::input_parameter< double >::type lphi(lphiSEXP);
+    rcpp_result_gen = Rcpp::wrap(compute_logk(mu, lmu, phi, lphi));
+    return rcpp_result_gen;
+END_RCPP
+}
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_flexcm_compute_logz", (DL_FUNC) &_flexcm_compute_logz, 2},
+    {"_flexcm_compute_logk", (DL_FUNC) &_flexcm_compute_logk, 4},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_flexcm(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}

src/normalizing_constants.cpp

@@ -0,0 +1,74 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+//' @title Compute COM-Poisson normalizing constante in log scale
+//' @description Computes the normalizing constant in the log scale
+//'   usign the LogSumExp trick to avoid numerical issues (See details).
+//' @details \code{logspace_add(a, b) = log(exp(a) + exp(b)) = b + log(1 +
+//'   exp(a - b))}, where \code{b > a}.
+//' @param loglambda A vector of logarithm of the \eqn{\lambda}
+//'   parameter.
+//' @param nu A vector of dispersion parameters \eqn{\nu}.
+//' @return The normalizing constant.
+//' @references Wikipedia, LogSumExp. \url{http://rstudio.com}.
+//' @author Eduardo Jr <[email protected]>
+//' @export
+// [[Rcpp::export]]
+NumericVector compute_logz(NumericVector loglambda,
+                           double nu) {
+  // Control loop
+  int maxiter = 1e4;
+  double logepsilon = log(1e-5);
+  // Output vector
+  int n = loglambda.size();
+  NumericVector out(n);
+  // Compute logz
+  for (int i = 0; i < n; ++i) {
+    double logz = 0;
+    for (int j = 1; j < maxiter; ++j) {
+      double logz_ = j * loglambda[i] - nu * R::lgammafn(j + 1);
+      logz = R::logspace_add(logz, logz_);
+      if (logz_ < logepsilon) break;
+    }
+    out[i] = logz;
+  }
+  return out;
+}
+
+//' @title Compute Double Poisson normalizing constante in log scale
+//' @description Computes the normalizing constant in the log scale
+//'   usign the LogSumExp trick to avoid numerical issues (See details).
+//' @details \code{logspace_add(a, b) = log(exp(a) + exp(b)) = b + log(1
+//'   + exp(a - b))}, where \code{b > a}.
+//' @param mu a vector of the \eqn{\mu} parameter.
+//' @param lmu logarithm of \eqn{\mu}.
+//' @param phi value of \eqn{\varphi} parameter.
+//' @param lphi logarithm of \eqn{\varphi}.
+//' @return The normalizing constant.
+//' @references Wikipedia, LogSumExp. \url{http://rstudio.com}.
+//' @author Eduardo Jr <[email protected]>
+//' @export
+// [[Rcpp::export]]
+NumericVector compute_logk(NumericVector mu,
+                           NumericVector lmu,
+                           double phi,
+                           double lphi) {
+  // Control loop
+  int maxiter = 1e4;
+  double logepsilon = log(1e-5);
+  // Output vector
+  int n = mu.size();
+  NumericVector out(n);
+  // Compute logz
+  for (int i = 0; i < n; ++i) {
+    double logk = -lphi/2 - mu[i]/phi;
+    for (int j = 1; j < maxiter; ++j) {
+      double logk_ = -0.5 * lphi - mu[i]/phi - j + j * log(j) -
+        R::lgammafn(j + 1) + (j/phi) * (1 + lmu[i] - log(j));
+      logk = R::logspace_add(logk, logk_);
+      if (logk_ < logepsilon) break;
+    }
+    out[i] = logk;
+  }
+  return out;
+}