Speed up null fitting for categorical X's via Fisher scoring + line search

R/dpseudohuber_median_d2x.R

@@ -0,0 +1,63 @@
+#' Calculate the second derivative of the pseudo-Huber smoothed median
+#'
+#' @param x A vector to calculate the derivative of the pseudo-Huber smoothed median for. 
+#' @param d Smoothing parameter, by default set to \code{0.1}. As \code{d} approaches 
+#' \code{0} the pseudo-Huber median function approaches the median and as \code{d} approaches 
+#' infinity this function approaches the mean. 
+#' @param na.rm Passed to \code{pseudohuber_median}, default is FALSE, if FALSE then when 
+#' \code{x} includes at least one NA value then NA is returned, if TRUE then when \code{x} 
+#' includes at least one NA value then that value is removed and the pseudo-Huber median 
+#' is computed without it.
+#'
+#' @return The second derivative of the calculated pseudo-Huber smoothed median over \code{x} with smoothing 
+#' parameter \code{d}. 
+#' 
+#' @note algebra derived by ChatGPT based on our supp
+#' 
+#' @examples
+#' dpseudohuber_median_dx(x = rnorm(10), d = 0.1)
+#'
+#' @aliases dpsuedohuber_median_dx
+#'
+#' @export
+#'
+dpseudohuber_median_d2x <- function(x,
+                                    d = 0.1,
+                                    na.rm = FALSE) {
+  # pseudo-Huber center
+  ps_center <- pseudohuber_median(x, d, na.rm = na.rm)
+
+  # basic quantities
+  y  <- (x - ps_center) / d
+  w  <- 1 / sqrt(1 + y^2)          # w_j
+  w3 <- w^3                        # a_j = w_j^3
+  S  <- sum(w3, na.rm = na.rm)     # S = sum_j w_j^3
+
+  # t_j = (x_j - c*)/d^2 * w_j^5
+  v  <- (x - ps_center) / d^2
+  w5 <- w^5
+  t  <- v * w5                     # t_j
+  T  <- sum(t, na.rm = na.rm)      # T = sum_j t_j
+
+  n <- length(x)
+  H <- matrix(NA_real_, n, n)
+
+  # Hessian entries:
+  # H_{pq} = -3/S * [ t_p (delta_{pq} - w3_q/S) - (t_q - w3_q*T/S) * w3_p/S ]
+  for (p in seq_len(n)) {
+    a_p <- w3[p]
+    t_p <- t[p]
+    for (q in seq_len(n)) {
+      a_q <- w3[q]
+      t_q <- t[q]
+      delta_pq <- as.numeric(p == q)
+
+      H[p, q] <- -(3 / S) * (
+        t_p * (delta_pq - a_q / S) -
+          (t_q - a_q * T / S) * a_p / S
+      )
+    }
+  }
+
+  H
+}

R/f_info.R

@@ -1,69 +1,69 @@
 #' @importFrom methods as
 # fisher information for multinomial model
-f_info <- function(Y,
-                   B_cup,
-                   B,
-                   X,
-                   X_cup,
-                   omit = NULL,
-                   compute_together = FALSE,
-                   return_separate = FALSE){
+f_info <- function(
+  Y,
+  B_cup,
+  B,
+  X,
+  X_cup,
+  omit = NULL,
+  compute_together = FALSE,
+  return_separate = FALSE
+) {
   n <- nrow(Y)
   J <- ncol(Y)
   p <- ncol(X)
 
-  if(is.null(omit)){
-  f_info <- matrix(0,ncol = p*J,nrow = p*J)
-  } else{
-    if(compute_together){
-    f_info <- matrix(0,ncol = p*(J - 1), nrow = p*(J - 1))
+  if (is.null(omit)) {
+    f_info <- matrix(0, ncol = p * J, nrow = p * J)
+  } else {
+    if (compute_together) {
+      f_info <- matrix(0, ncol = p * (J - 1), nrow = p * (J - 1))
     }
   }
-  z <- update_z(Matrix::Matrix(Y),Matrix::Matrix(X),B)
-
-
-
-  if(compute_together){ #older, slower way of computing info
-                        #used only to test equality with newer way
-    ident_mat <- Matrix::Diagonal(x = rep(1,J))
-    one_mat <- Matrix::Matrix(1,nrow = J,ncol=1)
-  for(i in 1:n){
-    which_rows <- 1:J + (i - 1)*J
-    mu_i <- exp(X_cup[which_rows,]%*% B_cup + z[i])
-    p_i <- mu_i/sum(mu_i)
-
-    M_iX_i <- (ident_mat - Matrix::tcrossprod(one_mat,p_i))%*%X_cup[which_rows,]
-
-
-    f_info <- f_info + Matrix::crossprod(M_iX_i,Matrix::Diagonal(x= as.numeric(mu_i)))%*%M_iX_i
-  }
-  } else{ #newer, faster way to compute info
-    mus <- exp(X_cup%*%B_cup + rep(z,each = J))
+  z <- update_z(Matrix::Matrix(Y), Matrix::Matrix(X), B)
+
+  if (compute_together) {
+    #older, slower way of computing info
+    #used only to test equality with newer way
+    ident_mat <- Matrix::Diagonal(x = rep(1, J))
+    one_mat <- Matrix::Matrix(1, nrow = J, ncol = 1)
+    for (i in 1:n) {
+      which_rows <- 1:J + (i - 1) * J
+      mu_i <- exp(X_cup[which_rows, ] %*% B_cup + z[i])
+      p_i <- mu_i / sum(mu_i)
+
+      M_iX_i <- (ident_mat - Matrix::tcrossprod(one_mat, p_i)) %*%
+        X_cup[which_rows, ]
+
+      f_info <- f_info +
+        Matrix::crossprod(M_iX_i, Matrix::Diagonal(x = as.numeric(mu_i))) %*%
+          M_iX_i
+    }
+  } else {
+    #newer, faster way to compute info
+    mus <- exp(X_cup %*% B_cup + rep(z, each = J))
     mu_diag <- Matrix::Diagonal(x = mus@x)
 
-    ps <- lapply(1:n,
-                 function(i){
-                   p_i <- mus@x[(i-1)*J +1:J];
-                 return(p_i/sum(p_i))})
-    ps_mat <- Matrix::sparseMatrix(i = 1:(n*J),
-                                   j = rep(1:n,each = J),
-                                   x = do.call(c,ps))
-    X_cup_p <- Matrix::crossprod(X_cup,ps_mat)
+    ps <- lapply(1:n, function(i) {
+      p_i <- mus@x[(i - 1) * J + 1:J]
+      return(p_i / sum(p_i))
+    })
+    ps_mat <- Matrix::sparseMatrix(
+      i = 1:(n * J),
+      j = rep(1:n, each = J),
+      x = do.call(c, ps)
+    )
+    X_cup_p <- Matrix::crossprod(X_cup, ps_mat)
     X_cup_p <- as.matrix(X_cup_p)
-    X_cup_p <- methods::as(X_cup_p,"sparseMatrix")
+    X_cup_p <- methods::as(X_cup_p, "sparseMatrix")
     sum_diag <- Matrix::Diagonal(x = rowSums(Y))
 
-    f_info <- Matrix::crossprod(X_cup,mu_diag)%*%X_cup -
-      X_cup_p%*%
-      Matrix::tcrossprod(sum_diag,X_cup_p)
-
-
-
-    }
-
-
+    f_info <- Matrix::crossprod(X_cup, mu_diag) %*%
+      X_cup -
+      X_cup_p %*%
+        Matrix::tcrossprod(sum_diag, X_cup_p)
+  }
 
   return(f_info)
-
-
 }