R: Как добавить джиттер только в единичных матрицах внутри функции?

Question

У меня есть следующая функция, которую мне нужно (m) применить к списку из более чем 1500 больших матриц (Z) и списка векторов (p) одинаковой длины. Тем не менее, я получаю ошибку, что некоторые матрицы являются единственными, как я уже разместил здесь. Вот моя функция:

kastner <- function(item, p) { print(item)
                               imp <- rowSums(Z[[item]])
                               exp <- colSums(Z[[item]])
                               x = p + imp
                               ac = p + imp - exp  
                               einsdurchx = 1/as.vector(x)    
                               einsdurchx[is.infinite(einsdurchx)] <- 0
                               A = Z[[item]] %*% diag(einsdurchx)
                               R = solve(diag(length(p))-A) %*% diag(p)    
                               C = ac * einsdurchx
                               R_bar = diag(as.vector(C)) %*% R
                               rR_bar = round(R_bar)
                               return(rR_bar)
}

и моя команда mapply, которая также печатает имена работающей матрицы:

KASTNER <- mapply(kastner, names(Z), p, SIMPLIFY = FALSE) 

Чтобы преодолеть проблему сингулярности, я хочу добавить небольшое количество jitter особые матрицы. Проблема начинается в строке 9 функции R = solve(diag(length(p))-A) %*% diag(p) как этот термин (diag(length(p))-A) становится единичным и не может быть solveд. Я попытался добавить джиттер для всех матриц Z в первой строке функции, используя: Z <- lapply(Z,function(x) jitter(x, factor = 0.0001, amount = NULL)), но это очень и очень мало и выдает все еще ошибки.

Поэтому моя идея состоит в том, чтобы проверить с if/else или что-то подобное if эта матрица diag(length(p))-A в единственном числе (возможно, используя собственные векторы для проверки коллинеарности) и добавить к этим матрицам дрожание, else (если нет) solve Команда должна выполняться как есть. Идеи, как реализовать это на функции? Спасибо

Вот некоторые примеры данных, хотя с сингулярностью проблем нет, так как я не смог перестроить эту ошибку для строки 9:

Z <- list("111.2012"= matrix(c(0,0,100,200,0,0,0,0,50,350,0,50,50,200,200,0),
                                 nrow = 4, ncol = 4, byrow = T),
          "112.2012"= matrix(c(10,90,0,30,10,90,0,10,200,50,10,350,150,100,200,10),
                                  nrow = 4, ncol = 4, byrow = T))
p <- list("111.2012"=c(200, 1000, 100, 10), "112.2012"=c(300, 900, 50, 100))

Изменить: небольшое количество дрожания не должно быть проблематичным в моих данных, так как у меня, вероятно, более 80% нулей в моих матрицах и чем большие значения. И меня интересуют только эти большие значения, но большое количество нулей, вероятно, является причиной сингулярности, но необходимо.

Answer 1

Поскольку вы не предоставили работающего примера, я не смог бы легко это проверить, поэтому бремя доказывания лежит на вас.:) В любом случае, это должно стать отправной точкой для дальнейшей работы. Комментарии в коде.

kastner <- function(item, p) { print(item)
  imp <- rowSums(Z[[item]])
  exp <- colSums(Z[[item]])
  x = p + imp
  ac = p + imp - exp  
  einsdurchx = 1/as.vector(x)    
  einsdurchx[is.infinite(einsdurchx)] <- 0

  # start a chunk that repeats until you get a valid result
  do.jitter <- TRUE # bureaucracy
  while (do.jitter == TRUE) {
    # run the code as usual
    A = Z[[item]] %*% diag(einsdurchx)
    # catch any possible errors, you can even catch "singularity" error here by
    # specifying error = function(e) e
    R <- tryCatch(solve(diag(length(p))-A) %*% diag(p), error = function(e) "jitterme")
    # if you were able to solve(), and the result is a matrix (carefuly if it's a vector!)...
    if (is.matrix(R)) {
      # ... turn the while loop off
      do.jitter <- FALSE
    } else {
      #... else apply some jitter and repeat by construcing A from a jittered Z[[item]]
      Z[[item]] <- jitter(Z[[item]])
    }
  }
  C = ac * einsdurchx
  R_bar = diag(as.vector(C)) %*% R
  rR_bar = round(R_bar)
  return(rR_bar)
}