Update documentation

thomvolker · Apr 5, 2024 · f0ee228 · f0ee228
1 parent eb765e6
commit f0ee228
Show file tree

Hide file tree

Showing 27 changed files with 895 additions and 218 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -10,16 +10,17 @@ S3method(print,kliep)
 S3method(print,lhss)
 S3method(print,naivedensityratio)
 S3method(print,naivesubspacedensityratio)
+S3method(print,spectral)
 S3method(print,summary.kliep)
 S3method(print,summary.lhss)
-S3method(print,summary.naivedensityratio)
-S3method(print,summary.naivesubspacedensityratio)
+S3method(print,summary.spectral)
 S3method(print,summary.ulsif)
 S3method(print,ulsif)
 S3method(summary,kliep)
 S3method(summary,lhss)
 S3method(summary,naivedensityratio)
 S3method(summary,naivesubspacedensityratio)
+S3method(summary,spectral)
 S3method(summary,ulsif)
 export(distance)
 export(kernel_gaussian)

diff --git a/R/permute.R b/R/permute.R
@@ -1,6 +1,6 @@
 #' Single permutation
 #' @rdname permute
-#' @param object Object of class \code{ulsif} or \code{kliep}
+#' @param object Density ratio object
 #' @param ... Additional arguments to pass through to specific permute functions.
 #' @return permutation statistic for a single permutation of the data
 #' @importFrom stats update
@@ -22,7 +22,7 @@ permute <- function(object, ...) {
 #' @importFrom stats predict
 #' @importFrom stats update
 
-permute.ulsif <- function(object, stacked, nnu, nde) {
+permute.ulsif <- function(object, stacked, nnu, nde, ...) {
   ind <- sample(rep(c(TRUE, FALSE), times = c(nnu, nde)))
   r <- stats::update(
     object = object,
@@ -37,7 +37,44 @@ permute.ulsif <- function(object, stacked, nnu, nde) {
   1/(2 * nnu) * sum(pred_nu) - 1/nde * sum(pred_de) + 1/2
 }
 
-permute.kliep <- function(object, stacked, nnu, nde) {
+#' Single permutation statistic of \code{kliep} object
+#' @rdname permute
+#' @param object \code{kliep} object
+#' @param stacked \code{matrix} with stacked numerator and denominator samples
+#' @param nnu Scalar with numerator sample size
+#' @param nde Scalar with denominator sample size
+#' @param min_pred Minimum value of the density ratio
+#' @return permutation statistic for a single permutation of the data
+#' @method permute kliep
+#' @importFrom stats predict
+#' @importFrom stats update
+
+permute.kliep <- function(object, stacked, nnu, nde, min_pred = sqrt(.Machine$double.eps), ...) {
+  ind <- sample(rep(c(TRUE, FALSE), times = c(nnu, nde)))
+  r <- stats::update(
+    object = object,
+    df_numerator = stacked[ind, ],
+    df_denominator = stacked[!ind, ],
+    progressbar = FALSE
+  )
+
+  pred_nu <- c(stats::predict(r, newdata = stacked[ind, , drop = FALSE]))
+
+  mean(log(pmax(min_pred, pred_nu)))
+}
+
+#' Single permutation statistic of \code{lhss} object
+#' @rdname permute
+#' @param object \code{lhss} object
+#' @param stacked \code{matrix} with stacked numerator and denominator samples
+#' @param nnu Scalar with numerator sample size
+#' @param nde Scalar with denominator sample size
+#' @return permutation statistic for a single permutation of the data
+#' @method permute lhss
+#' @importFrom stats predict
+#' @importFrom stats update
+
+permute.lhss <- function(object, stacked, nnu, nde, ...) {
   ind <- sample(rep(c(TRUE, FALSE), times = c(nnu, nde)))
   r <- stats::update(
     object = object,
@@ -47,11 +84,23 @@ permute.kliep <- function(object, stacked, nnu, nde) {
   )
 
   pred_nu <- c(stats::predict(r, newdata = stacked[ind, , drop = FALSE]))
+  pred_de <- c(stats::predict(r, newdata = stacked[!ind, , drop = FALSE]))
 
-  mean(log(pmax(sqrt(.Machine$double.eps), pred_nu)))
+  1/(2 * nnu) * sum(pred_nu) - 1/nde * sum(pred_de) + 1/2
 }
 
-permute.lhss <- function(object, stacked, nnu, nde) {
+#' Single permutation statistic of \code{spectral} object
+#' @rdname permute
+#' @param object \code{spectral} object
+#' @param stacked \code{matrix} with stacked numerator and denominator samples
+#' @param nnu Scalar with numerator sample size
+#' @param nde Scalar with denominator sample size
+#' @return permutation statistic for a single permutation of the data
+#' @method permute spectral
+#' @importFrom stats predict
+#' @importFrom stats update
+
+permute.spectral <- function(object, stacked, nnu, nde, ...) {
   ind <- sample(rep(c(TRUE, FALSE), times = c(nnu, nde)))
   r <- stats::update(
     object = object,
@@ -61,10 +110,24 @@ permute.lhss <- function(object, stacked, nnu, nde) {
   )
 
   pred_nu <- c(stats::predict(r, newdata = stacked[ind, , drop = FALSE]))
+  pred_de <- c(stats::predict(r, newdata = stacked[!ind, , drop = FALSE]))
 
-  mean(log(pmax(sqrt(.Machine$double.eps), pred_nu)))
+  1/(2 * nnu) * sum(pred_nu) - 1/nde * sum(pred_de) + 1/2
 }
 
+#' Single permutation statistic of \code{naivedensityratio} object
+#' @rdname permute
+#' @param object \code{naivedensityratio} object
+#' @param stacked \code{matrix} with stacked numerator and denominator samples
+#' @param nnu Scalar with numerator sample size
+#' @param nde Scalar with denominator sample size
+#' @param min_pred Minimum value of the predicted density ratio
+#' @param max_pred Maximum value of the predicted density ratio
+#' @return permutation statistic for a single permutation of the data
+#' @method permute naivedensityratio
+#' @importFrom stats predict
+#' @importFrom stats update
+
 permute.naivedensityratio <- function(object, stacked, nnu, nde, min_pred, max_pred) {
   ind <- sample(rep(c(TRUE, FALSE), times = c(nnu, nde)))
   r <- stats::update(
@@ -79,6 +142,19 @@ permute.naivedensityratio <- function(object, stacked, nnu, nde, min_pred, max_p
   (mean(pmin(max_pred, pmax(min_pred, pred_nu))) - mean(pmin(max_pred, pmax(min_pred, pred_de))))^2
 }
 
+#' Single permutation statistic of \code{naivesubspacedensityratio} object
+#' @rdname permute
+#' @param object \code{naivesubspacedensityratio} object
+#' @param stacked \code{matrix} with stacked numerator and denominator samples
+#' @param nnu Scalar with numerator sample size
+#' @param nde Scalar with denominator sample size
+#' @param min_pred Minimum value of the predicted density ratio
+#' @param max_pred Maximum value of the predicted density ratio
+#' @return permutation statistic for a single permutation of the data
+#' @method permute naivesubspacedensityratio
+#' @importFrom stats predict
+#' @importFrom stats update
+
 permute.naivesubspacedensityratio <- function(object, stacked, nnu, nde, min_pred, max_pred) {
   ind <- sample(rep(c(TRUE, FALSE), times = c(nnu, nde)))
   r <- stats::update(

diff --git a/R/predict.R b/R/predict.R
@@ -141,8 +141,6 @@ predict.lhss <- function(object, newdata = NULL, sigma = c("sigmaopt", "all"), l
 #' @param newdata Optional \code{matrix} new data set to compute the density
 #' @param sigma A scalar with the Gaussian kernel width
 #' @param J integer indicating the dimension of the eigenvector expansion
-#' @param tol A scalar indicating the smallest eligible density ratio value
-#' (used to censor negative predicted density ratio values).
 #' @param ... Additional arguments to be passed to the function
 #'
 #' @return An array with predicted density ratio values from possibly new data,
@@ -162,7 +160,7 @@ predict.lhss <- function(object, newdata = NULL, sigma = c("sigmaopt", "all"), l
 #' predict(fit1, newdata = rbind(x, y), sigma = 2, J = 10)
 
 predict.spectral <- function(object, newdata = NULL, sigma = c("sigmaopt", "all"),
-                             J = c("Jopt", "all"), tol = 1e-6, ...) {
+                             J = c("Jopt", "all"), ...) {
 
   newsigma <- check.sigma.predict(object, sigma)
   newJ     <- check.J.predict(object, J)
@@ -173,9 +171,11 @@ predict.spectral <- function(object, newdata = NULL, sigma = c("sigmaopt", "all"
 
   for (i in 1:length(newsigma)) {
     K <- distance(newdata, object$centers, FALSE) |> kernel_gaussian(newsigma[i])
-    phihatpred <- K %*% alpha_eigen$Evecs[,,i] %*% diag(sqrt(nrow(object$centers))/alpha_eigen$Evals[,i])
+    D <- diag(length(alpha_eigen$Evals[,i]))
+    diag(D) <- sqrt(nrow(object$centers))/alpha_eigen$Evals[,i]
+    phihatpred <- K %*% alpha_eigen$Evecs[,,i] %*% D
     for (j in 1:length(newJ)) {
-      dratio[ , j, i] <- pmax(tol, phihatpred[,seq_len(newJ[j])] %*% alpha_eigen$alpha[seq_len(newJ[j]),i])
+      dratio[ , j, i] <- phihatpred[,seq_len(newJ[j]), drop = FALSE] %*% alpha_eigen$alpha[seq_len(newJ[j]),i, drop = FALSE]
     }
   }
   dratio