R/get_dists.R

In devoges/fstat: fstat

#' Title
#'
#' @param Y n by 1 vector
#' @param X n by p matrix
#' @param Zs list of n by i_k matrices. NAMED
#' @param Ss list of n by n matrices, NAMED
#' @param prefix what to append to beginning of messages
#' @param print whether to print information
#' @param vcs_HA vector of variance component estimates under alternative
#' @param test which element of Zs or Ss to test
#' @param samples how many samples
#'
#' @return
#' @export
#'
#' @examples
get_dists <- function(Y,
                      X,
                      Zs = NULL,
                      Ss = NULL,
                      vcs_HA,
                      test,
                      prefix = NULL,
                      print = TRUE,
                      samples = 10000) {
  ##-- get basic things
  n <- X %>% nrow
  p <- X %>% ncol
  q <- p
  Xt <- t(X)
  effects <- Zs %>% names %>% c(Ss %>% names)
  nuis <- effects %>% remove_these(test)

  ##-- Fix gamma at value under H0
  gamma <- vcs_HA[nuis]

  ##-- Get anhilator matrix of X
  A <- X %mtimes% minv(Xt %mtimes% X) %mtimes% Xt

  ##-- Turn Ss into Zs
  Zs_from_Ss <- nlapply(Ss %>% names, function(x) {
    t(chol(Ss[[x]]))
  })
  Zs %<>% c(Zs_from_Ss)

  ##-- Save Zts
  Zts <- nlapply(effects, function(n) {
    t(Zs[[n]])
  })

  ##-- Get Sigma^0.5, they are all diagonal
  Sigma0.5s <- nlapply(effects, function(n) {
    Zs[[n]] %>% ncol %>% diag
  })

  ##-- Let C_n = Z_n \times Sigma^0.5_k
  message(prefix %% ": calculating Cs...")
  tic()
  Cs <- nlapply(effects, function(n) {
    Zs[[n]] %mtimes% Sigma0.5s[[n]]
  })
  Cts <- nlapply(effects, function(n) {
    t(Cs[[n]])
  })
  done <- toc(quiet = TRUE)
  time <- done$toc - done$tic
  message(prefix %% ": calculating Cs took " %% time %% " seconds")

  message(prefix %% ": calculating Ks...")
  tic()
  ##-- Let K = C_k \times C_k'
  Ks <- nlapply(effects, function(n) {
    Cs[[n]] %mtimes% Cts[[n]]
  })
  done <- toc(quiet = TRUE)
  time <- done$toc - done$tic
  message(prefix %% ": calculating Ks took " %% time %% " seconds")

  ##-- Function to calculate V_0(\gamma)^-1
  v0inv_gamma <- function(gamma) {
    message(prefix %% ": calculating v0inv_gamma...")
    tic()
    mats <- nlapply(nuis, function(n) {
      gamma[n] * Ks[[n]]
    })
    out <- minv(diag(n) + Reduce(`+`, mats))
    done <- toc(quiet = TRUE)
    time <- done$toc - done$tic
    message(prefix %% ": calculating v0inv_gamma took " %% time %% " seconds")
    return(out)
  }

  ##-- Function to calculate PV_0(\gamma)
  Pv0_gamma <- function(gamma) {
    message(prefix %% ": calculating Pv0_gamma...")
    tic()
    vinv <- v0inv_gamma(gamma)
    vinvX <- vinv %mtimes% X
    middle <- minv(Xt %mtimes% vinv %mtimes% X)
    second <- vinvX %mtimes% middle %mtimes% t(vinvX)
    out <- vinv - second
    done <- toc(quiet = TRUE)
    time <- done$toc - done$tic
    message(prefix %% ": calculating Pv0_gamma took " %% time %% " seconds")
    return(out)
  }

  ##-- Function to calculate \rho(\gamma)
  rho_gamma <- function(gamma) {
    message(prefix %% ": calculating rho...")
    tic()
    Pv0 <- Pv0_gamma(gamma)
    out <- meig(Cts[[test]] %mtimes% Pv0 %mtimes% Cs[[test]])
    done <- toc(quiet = TRUE)
    time <- done$toc - done$tic
    message(prefix %% ": calculating rho took " %% time %% " seconds")
    return(out)
  }

  ##-- calculate rho
  if ( length(gamma) > 0 ) {
    rho <- rho_gamma(gamma)
  } else {
    ## no nuis parameters
    message(prefix %% ": calculating rho...no nuissance parameters")
    tic()
    v0 <- diag(n)
    v0inv <- minv(v0)
    P <- minv(v0) - X %mtimes% minv(Xt %mtimes% v0inv %mtimes% X) %mtimes% Xt %mtimes% v0inv
    rho <- meig(Cs[[test]] %mtimes% P %mtimes% Cts[[test]])
    done <- toc(quiet = TRUE)
    time <- done$toc - done$tic
    message(prefix %% ": calculating rho took " %% time %% " seconds")
  }

  ##-- remove rho that are basically 0
  if ( length(rho) > (n - p) ) {
    rho <- rho[1:(n - p)]
  }

  ##-- get omega eigen values
  message(prefix %% ": calculating omega eigen values...")
  tic()
  P1 <- diag(n) - A
  ws <- nlapply(nuis, function(n) {
    meig(Cts[[n]] %mtimes% P1 %mtimes% Cs[[n]])
  })
  done <- toc(quiet = TRUE)
  time <- done$toc - done$tic
  message(prefix %% ": calculating omega eigen values took " %% time %% " seconds")

  checks <- quantile(1:samples, probs = seq(0.1, 0.9, 0.1)) %>% round
  one_sample <- function(x) {
    if ( any(x == checks) & print ) {
      e <- which(x == checks) %>% names % % "of samples done"
      message(prefix %% ": " %% e)
    }

    u2 <- rnorm(n - p)^2
    v2 <- rchisq(1, df = p - q)

    to_max_f <- function(lambda) {
      fn(u2 = u2,
         lambda = lambda,
         rho = rho,
         w = ws,
         gamma = gamma)
    }

    to_max_h <- function(lambda) {
      hn(u2 = u2,
           lambda = lambda,
           rho = rho)
    }

    for_rlrt <- optimise(f = to_max_h,
                            interval = c(0, 100),
                            maximum = TRUE,
                            tol = .Machine$double.eps)

    rlrt_lambda <- for_rlrt$maximum
    rlrt <- for_rlrt$objective

    f_lambda <- optimise(f = to_max_f,
                       interval = c(0, 1000),
                       maximum = TRUE,
                       tol = .Machine$double.eps)$maximum
    f <- fstat(u2 = u2,
               v2 = v2,
               lrho = f_lambda * rho)

    tibble(Type = c("f", "rlrt"),
           Sample = c(f, rlrt),
           Lambda = c(f_lambda, rlrt_lambda))
  }

  message(prefix %% ": finally calculating samples...")
  tic()
  info <- lapply(1:samples, one_sample) %>% bind_rows
  done <- toc(quiet = TRUE)
  time <- done$toc - done$tic
  message(prefix %% ": calculating samples took " %% time %% " seconds")


  out <- list()
  out$info <- info
  out$rho <- rho
  return(out)
}