Add C++ log-likelihood function to speed up optimization in dev_beta_binom()

DESCRIPTION

@@ -26,6 +26,7 @@ Depends:
 Imports:
     chk,
     lifecycle,
+    Rcpp (>= 1.0.12),
     stats
 Suggests: 
     aods3,
@@ -43,6 +44,7 @@ Suggests:
     tidyr,
     viridis,
     withr
+LinkingTo: Rcpp
 Config/testthat/edition: 3
 Encoding: UTF-8
 Language: en-US

NAMESPACE

@@ -150,8 +150,10 @@ export(xtr_sd)
 export(zeros)
 export(zscore)
 import(chk)
+importFrom(Rcpp,sourceCpp)
 importFrom(stats,dbinom)
 importFrom(stats,dlnorm)
 importFrom(stats,dnbinom)
 importFrom(stats,dnorm)
 importFrom(stats,dpois)
+useDynLib(extras)

R/RcppExports.R

@@ -0,0 +1,7 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+log_lik_beta_binom2 <- function(x, size, prob, theta) {
+    .Call('_extras_log_lik_beta_binom2', PACKAGE = 'extras', x, size, prob, theta)
+}
+

R/dev.R

@@ -13,7 +13,7 @@
 #' dev_beta_binom(c(0, 1, 2), 1, 0.5, 0)
 dev_beta_binom <- function(x, size = 1, prob = 0.5, theta = 0, res = FALSE) {
   opt_beta_binom <- function(prob, x, size = size, theta = theta) {
-    -log_lik_beta_binom(x = x, size = size, prob = prob, theta = theta)
+    -log_lik_beta_binom2(x = x, size = size, prob = prob, theta = theta)
   }
   if (length(size) == 1) {
     size <- rep(size, length(x))

R/namespace.R

@@ -1,3 +1,5 @@
 #' @import chk
 #' @importFrom stats dbinom dlnorm dnorm dpois dnbinom
+#' @useDynLib extras
+#' @importFrom Rcpp sourceCpp
 NULL

src/RcppExports.cpp

@@ -0,0 +1,36 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+#ifdef RCPP_USE_GLOBAL_ROSTREAM
+Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
+#endif
+
+// log_lik_beta_binom2
+double log_lik_beta_binom2(int x, int size, double prob, double theta);
+RcppExport SEXP _extras_log_lik_beta_binom2(SEXP xSEXP, SEXP sizeSEXP, SEXP probSEXP, SEXP thetaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int >::type x(xSEXP);
+    Rcpp::traits::input_parameter< int >::type size(sizeSEXP);
+    Rcpp::traits::input_parameter< double >::type prob(probSEXP);
+    Rcpp::traits::input_parameter< double >::type theta(thetaSEXP);
+    rcpp_result_gen = Rcpp::wrap(log_lik_beta_binom2(x, size, prob, theta));
+    return rcpp_result_gen;
+END_RCPP
+}
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_extras_log_lik_beta_binom2", (DL_FUNC) &_extras_log_lik_beta_binom2, 4},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_extras(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}

src/log_lik_beta_binom.cpp

@@ -0,0 +1,34 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// @noRd
+// [[Rcpp::export]]
+double log_lik_beta_binom2(int x, int size, double prob, double theta) {
+  double alpha = prob * 2 * (1 / theta);
+  double beta = (1 - prob) * 2 * (1 / theta);
+  double lbeta_binom = R::lgammafn(size + 1) - R::lgammafn(x + 1) - R::lgammafn(size - x + 1) +
+    R::lgammafn(x + alpha) + R::lgammafn(size - x + beta) - R::lgammafn(size + alpha + beta) +
+    R::lgammafn(alpha + beta) - R::lgammafn(alpha) - R::lgammafn(beta);
+
+  bool use_binom = theta == 0;
+
+  if (use_binom) {
+    return(R::dbinom(x, size, prob, 1));
+  }
+  if ((((x == 0) & (prob == 0)) | ((x == size) & (prob == 1)))) {
+    return(0.0);
+  }
+  if ((x != 0) & (prob == 0)) {
+    return(R_NegInf);
+  }
+  if ((x != size) & (prob == 1)) {
+    return(R_NegInf);
+  }
+  if ((x > size)) {
+    return(R_NegInf);
+  }
+  if ((theta < 0)) {
+    return(R_NaN);
+  }
+  return(lbeta_binom);
+}