retel() wip

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
 export(etel)
+export(retel)
 importFrom(Matrix,rankMatrix)
 importFrom(nloptr,nloptr)

R/etel.R

@@ -2,41 +2,50 @@
 #'
 #' Computes exponentially tilted empirical likelihood.
 #'
-#' @param g A numeric matrix, or an object that can be coerced to a numeric
-#'   matrix. Each row corresponds to an observation of an estimating function.
-#'   The number of rows must be greater than the number of columns.
+#' @param fn An estimating function that takes the data `x` and parameter value
+#'   `par` as its arguments, returning a numeric matrix. Each row is the return
+#'   value from the corresponding row in `x`.
+#' @param x A numeric matrix, or an object that can be coerced to a numeric
+#'   matrix. Each row corresponds to an observation. The number of rows must be
+#'   greater than the number of columns.
+#' @param par A numeric vector of parameter values to be tested. The length of
+#'   the vector must be the same as the number of columns in `x`.
 #' @param opts A list with optimization options for [nloptr()].
 #' @return A single numeric of the log-likelihood ratio.
 #' @references Schennach, SM (2005).
 #'   "Bayesian exponentially tilted empirical likelihood"
 #'   \emph{Biometrika}, 92, 31--46.
 #' @examples
 #' set.seed(63456)
+#' f <- function(x, par) {
+#'   x - par
+#' }
 #' x <- rnorm(100)
 #' par <- 0
-#' g <- x - par
-#' etel(g)
+#' etel(f, x, par)
 #' @export
-etel <- function(g, opts) {
-  mm <- as.matrix(g, rownames.force = TRUE)
-  n <- nrow(mm)
-  p <- ncol(mm)
+etel <- function(fn, x, par, opts) {
+  x <- validate_x(x)
+  par <- validate_par(par)
+  g <- fn(x, par)
+  g <- as.matrix(g, rownames.force = TRUE)
+  n <- nrow(g)
+  p <- ncol(g)
   stopifnot(
     "`g` must have at least two observations." = (n >= 2L),
     "`g` must be a finite numeric matrix." =
-      (isTRUE(is.numeric(mm) && all(is.finite(mm)))),
-    "`g` must have full column rank." = (isTRUE(n > p && rankMatrix(mm) == p))
+      (isTRUE(is.numeric(g) && all(is.finite(g)))),
+    "`g` must have full column rank." = (isTRUE(n > p && rankMatrix(g) == p))
   )
   if (missing(opts)) {
     opts <- list("algorithm" = "NLOPT_LD_LBFGS", "xtol_rel" = 1e-06)
   }
   optim <- nloptr(
-    x0 = rep(0, NCOL(g)), eval_f = eval_d_fun, eval_grad_f = eval_gr_d_fun,
-    opts = opts, g = mm, n = n
+    x0 = rep(0, ncol(g)), eval_f = eval_d_fn, eval_grad_f = eval_gr_d_fn,
+    opts = opts, g = g, n = n
   )
   lambda <- optim$solution
-  out <-
-    -n * log(eval_d_fun(lambda, mm, n)) + as.numeric(lambda %*% colSums(mm))
+  out <- as.numeric(lambda %*% colSums(g)) - n * log(eval_d_fn(lambda, g, n))
   attributes(out) <- list(optim = optim)
   out
 }

R/helpers.R

@@ -7,10 +7,41 @@
 #' @param n A single integer.
 #' @return A single numeric.
 #' @noRd
-eval_d_fun <- function(l, g, n) {
+eval_d_fn <- function(l, g, n) {
   sum(exp(g %*% l)) / n
 }
 
+#' Exponentially tilted empirical likelihood objective function
+#'
+#' Computes the exponentially tilted empirical likelihood objective function.
+#'
+#' @param l A numeric vector.
+#' @param mu A numeric matrix.
+#' @param Sigma A numeric matrix.
+#' @param tau A single numeric.
+#' @param n A single integer.
+#' @return A single numeric.
+#' @noRd
+eval_p_fn <- function(l, mu, Sigma, tau, n) {
+  tau * exp(as.numeric(l %*% mu) + colSums(mu * (Sigma %*% mu)) / 2) / n
+}
+
+#' Exponentially tilted empirical likelihood objective function
+#'
+#' Computes the exponentially tilted empirical likelihood objective function.
+#'
+#' @param l A numeric vector.
+#' @param g A numeric matrix.
+#' @param mu A numeric matrix.
+#' @param Sigma A numeric matrix.
+#' @param tau A single numeric.
+#' @param n A single integer.
+#' @return A single numeric.
+#' @noRd
+eval_obj_fn <- function(l, g, mu, Sigma, tau, n) {
+  eval_d_fn(l, g, n) + eval_p_fn(l, mu, Sigma, tau, n)
+}
+
 #' Gradient of exponentially tilted empirical likelihood objective function
 #'
 #' Computes the gradient of the exponentially tilted empirical likelihood
@@ -19,8 +50,77 @@ eval_d_fun <- function(l, g, n) {
 #' @param l A numeric vector.
 #' @param g A numeric matrix.
 #' @param n A single integer.
-#' @return A single numeric.
+#' @return A numeric vector.
+#' @noRd
+eval_gr_d_fn <- function(l, g, n) {
+  colSums(as.numeric(exp(g %*% l)) * g) / n
+}
+
+#' Gradient of exponentially tilted empirical likelihood objective function
+#'
+#' Computes the gradient of the exponentially tilted empirical likelihood
+#' objective function.
+#'
+#' @param l A numeric vector.
+#' @param mu A numeric matrix.
+#' @param Sigma A numeric matrix.
+#' @param tau A single numeric.
+#' @param n A single integer.
+#' @return A numeric vector.
+#' @noRd
+eval_gr_p_fn <- function(l, mu, Sigma, tau, n) {
+  eval_p_fn(l, mu, Sigma, tau, n) * (mu + as.numeric(Sigma %*% l))
+}
+
+#' Gradient of exponentially tilted empirical likelihood objective function
+#'
+#' Computes the gradient of the exponentially tilted empirical likelihood
+#' objective function.
+#'
+#' @param l A numeric vector.
+#' @param g A numeric matrix.
+#' @param mu A numeric matrix.
+#' @param Sigma A numeric matrix.
+#' @param tau A single numeric.
+#' @param n A single integer.
+#' @return A numeric vector.
+#' @noRd
+eval_gr_obj_fn <- function(l, g, mu, Sigma, tau, n) {
+  eval_gr_d_fn(l, g, n) + eval_gr_p_fn(l, mu, Sigma, tau, n)
+}
+
+#' Validate `x`
+#'
+#' Validate `x` in [etel()] and its extensions.
+#'
+#' @param x A numeric matrix, or an object that can be coerced to a numeric
+#'   matrix.
+#' @return A numeric matrix.
+#' @noRd
+validate_x <- function(x) {
+  x <- as.matrix(x, rownames.force = TRUE)
+  stopifnot(
+    "`x` must have at least two observations." = (nrow(x) >= 2L),
+    "`x` must must have larger number of rows than columns." =
+      nrow(x) > ncol(x),
+    "`x` must be a finite numeric matrix." =
+      isTRUE(is.numeric(x) && all(is.finite(x))),
+    "`x` must have full column rank." = rankMatrix(x) == ncol(x)
+  )
+  x
+}
+
+#' Validate `par`
+#'
+#' Validate `par` in [etel()] and its extensions.
+#'
+#' @param par A numeric vector.
+#' @return A numeric vector.
 #' @noRd
-eval_gr_d_fun <- function(l, g, n) {
-  colSums(as.vector(exp(g %*% l)) * g) / n
+validate_par <- function(par) {
+  stopifnot(
+    "`par` must be a finite numeric vector." =
+      (isTRUE(is.numeric(par) && all(is.finite(par))))
+  )
+  par
 }

R/retel-package.R

@@ -1,3 +1,4 @@
+#' @aliases retel-package
 #' @references Kim E, MacEachern SN, Peruggia M (2023).
 #'   "Regularized Exponentially Tilted Empirical Likelihood for Bayesian
 #'   Inference."

R/retel.R

@@ -0,0 +1,61 @@
+#' Regularized exponentially tilted empirical likelihood
+#'
+#' Computes regularized exponentially tilted empirical likelihood.
+#'
+#' @param fn An estimating function that takes the data `x` and parameter value
+#'   `par` as its arguments, returning a numeric matrix. Each row is the return
+#'   value from the corresponding row in `x`.
+#' @param x A numeric matrix, or an object that can be coerced to a numeric
+#'   matrix. Each row corresponds to an observation. The number of rows must be
+#'   greater than the number of columns.
+#' @param par A numeric vector of parameter values to be tested. The length of
+#'   the vector must be the same as the number of columns in `x`.
+#' @param mu A numeric matrix.
+#' @param Sigma A numeric matrix.
+#' @param tau A single numeric.
+#' @param opts A list with optimization options for [nloptr()].
+#' @return A single numeric of the log-likelihood ratio.
+#' @references Kim E, MacEachern SN, Peruggia M (2023).
+#'   "Regularized Exponentially Tilted Empirical Likelihood for Bayesian
+#'   Inference."
+#'   \doi{10.48550/arXiv.2312.17015}.
+#' @examples
+#' set.seed(63456)
+#' f <- function(x, par) {
+#'   x - par
+#' }
+#' x <- rnorm(100)
+#' par <- 0
+#' mu <- 0
+#' Sigma <- matrix(rnorm(1), nrow = 1)
+#' tau <- 1
+#' retel(f, x, par, mu, Sigma, tau)
+#' @export
+retel <- function(fn, x, par, mu, Sigma, tau, opts) {
+  x <- validate_x(x)
+  par <- validate_par(par)
+  g <- fn(x, par)
+  g <- as.matrix(g, rownames.force = TRUE)
+  n <- nrow(g)
+  p <- ncol(g)
+  stopifnot(
+    "`g` must have at least two observations." = (n >= 2L),
+    "`g` must be a finite numeric matrix." =
+      (isTRUE(is.numeric(g) && all(is.finite(g)))),
+    "`g` must have full column rank." = (isTRUE(n > p && rankMatrix(g) == p))
+  )
+  if (missing(opts)) {
+    opts <- list("algorithm" = "NLOPT_LD_LBFGS", "xtol_rel" = 1e-06)
+  }
+  optim <- nloptr(
+    x0 = rep(0, ncol(g)), eval_f = eval_obj_fn, eval_grad_f = eval_gr_obj_fn,
+    opts = opts, g = g, mu = mu, Sigma = Sigma, tau = tau, n = n
+  )
+  lambda <- optim$solution
+  out <- as.numeric(lambda %*% colSums(g)) + as.numeric(lambda %*% mu) +
+    colSums(mu * (Sigma %*% mu)) / 2 -
+    (n + 1) * log(n) + (n + 1) * log(n + tau) -
+    (n + 1) * log(eval_obj_fn(lambda, g, mu, Sigma, tau, n))
+  attributes(out) <- list(optim = optim)
+  out
+}