R/cond_BIC.R

In haodongucsb/toyexample:

#' for nei
#' @import quadprog
#' @import gss
#' @export
#' @param Dmat Matrix
#' @param dvec vector
#' @param Amat vector
#' @param bvec vector
My_solve.QP4 <- function(Dmat, dvec, Amat, bvec) {
  solution <- tryCatch(solve.QP(Dmat, dvec, Amat, bvec)$solution, error = function(x) NA)
  if (is.na(solution[1])) {
    M <- solve(Dmat)
    Dmat <- t(M) %*% M
    sc <- norm(Dmat, "2")
    solution <- tryCatch(solve.QP(Dmat = Dmat / sc, dvec = dvec / sc, Amat = Amat, bvec = bvec, meq = 0, factorized = FALSE)$solution, error = function(x) NA)
    if (is.na(solution[1])) {
      Dmat <- diag(diag(Dmat))
      solution <- solve.QP(Dmat, dvec, Amat, bvec)$solution
    }
  }
  return(solution)
}
cond_BIC <- function(y, formula, response, type = NULL, ydomain = as.list(NULL), M0, M_list, maxiteration, tolerance, id.basis = NULL, theta1 = NULL, theta2, w2 = NULL, n, dimen, p_f = 0) {
  p <- dimen
  loop <- 1
  ltheta2 <- length(theta2)
  #c_matrix <- matrix(0, c_length, maxiteration)
  Theta2 <- matrix(0, ltheta2, maxiteration)
  Theta2[, loop] <- theta2
  Theta1 <- matrix(0, p, maxiteration)
  l <- rep(0, maxiteration)
  l[loop] <- 0
  diff2 <- 1
  d <- rep(0, maxiteration)
  d[1] <- diff2
  zerodiff <- 1
  check_na <- 0
  while (diff2 >= tolerance & loop < maxiteration & (zerodiff) >= 1) {
    loop <- loop + 1
    densityfit <- sscden_selection(formula = formula, response = response, data = y, ydomain = ydomain, id.basis = id.basis, theta2 = theta2, w2 = w2, seed = 5732, p = p, p_f = p_f, type = type)
    M <- cond_select_M(densityfit, w2, M0, M_list, k = 5, n, p = p, p_f = p_f)[[1]]
    theta1 <- densityfit$theta1
    theta2.1 <- densityfit$theta2
    R <- densityfit$R
    ltheta2 <- length(theta2.1)
    c1 <- densityfit$c
    d1 <- densityfit$d
    U <- densityfit$r
    U2 <- NULL
    for (i in (p + 1 - p_f):(2 * p - 1 - p_f)) {
      U2 <- cbind(U2, U[, , i])
    }
    u <- densityfit$rbasis
    u2 <- NULL
    for (i in (p + 1 - p_f):(2 * p - 1 - p_f)) {
      u2 <- cbind(u2, u[, , i])
    }
    g <- densityfit$g
    t2 <- rep(0, length(g))
    for (i in 1:length(g)) {
      t2[i] <- exp(-g[i]) / sum(exp(-g))
    }
    B2 <- densityfit$int.r[, (p + 1 - p_f):(2 * p - 1 - p_f)]
    la1 <- 10^densityfit$lambda / 2
    Hpart1 <- U2 %*% kronecker(diag(ltheta2), c1)
    Gmatrix <- -t(Hpart1) %*% t2 + t(B2) %*% c1 + la1 * kronecker(diag(ltheta2), t(c1)) %*% t(u2) %*% c1
    Hpart2 <- t(t2) %*% U2 %*% kronecker(diag(ltheta2), c1)
    Hpart3 <- sqrt(t2) * Hpart1
    Hmatrix <- t(Hpart3) %*% Hpart3 - t(Hpart2) %*% Hpart2
    if (is.positive.definite(Hmatrix)) {
      Hmatrix <- Hmatrix
    } else {
      Hmatrix <- Hmatrix + diag(rep(1e-6, dim(Hmatrix)[1]))
    }
    dvec <- Hmatrix %*% theta2.1 - Gmatrix
    Dmat <- Hmatrix
    Amat <- rbind(diag(ltheta2), -w2)
    bvec <- c(rep(0, ltheta2), -M)
    theta2 <- My_solve.QP(Dmat, dvec, t(Amat), bvec)
    theta2[theta2 < 1e-10] <- 0
    l[loop] <- la1
    Theta2[, loop] <- theta2
    if (loop <= 3) {
      zerodiff <- 1
    } else {
      zerodiff <- sum(theta2 > 0) - sum(Theta2[, loop - 3] > 0)
    }
    diff2 <- sqrt(sum((Theta2[, loop] - Theta2[, loop - 1])^2)) / (sqrt(sum((Theta2[, loop - 1])^2)) + 1e-6)
    #diff_c <- sqrt(sum((c_matrix[, loop] - c_matrix[, loop - 1])^2)) / (sqrt(sum((c_matrix[, loop - 1])^2)) + 1e-6)
    d[loop] <- diff2
  }

  thre_list <- c(1e-2, 1e-3, 1e-4)
  count <- 0
  thre_mat <- rep(0, length(thre_list))
  for (thre in thre_list) {
    count <- count + 1
    theta20 <- theta2
    theta20[theta20 < thre] <- 0
    R0 <- matrix(0, dim(densityfit$r)[1], dim(densityfit$r)[2])
    for (i in 1:(p - p_f)) {
      R0 <- R0 + 10^densityfit$theta1[i] * densityfit$r[, , i]
    }
    for (i in (p + 1 - p_f):(2 * p - 1 - p_f)) {
      R0 <- R0 + theta20[i - p + p_f] * densityfit$r[, , i] / w2[i - p + p_f]
    }
    cv_g <- densityfit$s %*% d1 + R0 %*% c1
    score <- (sum(exp(-cv_g)))
    thre_mat[count] <- log(score / n) + densityfit$loss_int
  }
  min_pos <- which.min(thre_mat)
  bound <- thre_list[min_pos]
  theta2[theta2 < bound] <- 0
  list(loop = loop, theta2 = theta2, densityfit = densityfit)
}


cond_select_M <- function(densityfit, w2, M0, M_list, k = 1, n, p, p_f = 0) {
  M_list1 <- M0 * M_list
  cv_mat <- rep(0, length(M_list1))
  count <- 0
  densityfit <- densityfit
  theta1.1 <- densityfit$theta1
  theta2.1 <- densityfit$theta2
  ltheta2 <- length(theta2.1)
  c1 <- densityfit$c
  d1 <- densityfit$d
  R <- densityfit$R
  ltheta2 <- length(theta2.1)
  U <- densityfit$r
  int.r <- densityfit$int.r
  U2 <- NULL
  for (i in (p + 1 - p_f):(2 * p - 1 - p_f)) {
    U2 <- cbind(U2, U[, , i])
  }
  u <- densityfit$rbasis
  u2 <- NULL
  for (i in (p + 1 - p_f):(2 * p - 1 - p_f)) {
    u2 <- cbind(u2, u[, , i])
  }
  g <- densityfit$g
  t2 <- rep(0, length(g))
  for (i in 1:length(g)) {
    t2[i] <- exp(-g[i]) / sum(exp(-g))
  }
  B2 <- densityfit$int.r[, (p + 1 - p_f):(ltheta2 + p - p_f)]
  la1 <- 10^densityfit$lambda / 2
  Hpart1 <- U2 %*% kronecker(diag(ltheta2), c1)
  Gmatrix <- -t(Hpart1) %*% t2 + t(B2) %*% c1 + la1 * kronecker(diag(ltheta2), t(c1)) %*% t(u2) %*% c1
  Hpart2 <- t(t2) %*% U2 %*% kronecker(diag(ltheta2), c1)
  Hpart3 <- sqrt(t2) * Hpart1
  Hmatrix <- t(Hpart3) %*% Hpart3 - t(Hpart2) %*% Hpart2
  if (is.positive.definite(Hmatrix)) {
    Hmatrix <- Hmatrix
  } else {
    Hmatrix <- Hmatrix + diag(rep(1e-6, dim(Hmatrix)[1]))
  }
  dvec <- Hmatrix %*% theta2.1 - Gmatrix
  Dmat <- Hmatrix
  Amat <- rbind(diag(ltheta2), -w2)
  for (M in M_list1) {
    count <- count + 1
    bvec <- c(rep(0, ltheta2), -M)
    theta2 <- My_solve.QP(Dmat, dvec, t(Amat), bvec)
    theta2[theta2 < 1e-10] <- 0
    int.r.wk <- 0
    for (i in 1:(p - p_f)) {
      int.r.wk <- int.r.wk + 10^theta1.1[i] * int.r[, i]
    }
    for (i in (p + 1 - p_f):(2 * p - 1 - p_f)) {
      int.r.wk <- int.r.wk + theta2[i - p + p_f] * int.r[, i] / w2[i - p + p_f]
    }
    r <- densityfit$r
    R0 <- densityfit$R1
    for (i in (p + 1 - p_f):(2 * p - 1 - p_f)) {
      R0 <- R0 + theta2[i - p + p_f] * r[, , i] / w2[i - p + p_f]
    }
    cv_g <- densityfit$s %*% d1 + R0 %*% c1
    score <- (sum(exp(-cv_g)))
    cv_mat[count] <- log(score / n) + (dot(int.r.wk, c1) + dot(densityfit$int.s, d1)) / n + log(n * sum(theta2 > 1e-10))
  }
  min_error_pos <- which.min(cv_mat)
  M_opt <- M_list1[min_error_pos]
  list(M_opt = M_opt)
}