1

我有以下功能，我需要 (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)因为这个 term( diag(length(p))-A) 变得单数并且不能是solved。我尝试在函数的第一行使用: 将抖动添加到所有 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))

编辑：少量 o 抖动在我的数据中不应该有问题，因为我的矩阵中可能有超过 80% 的零而不是大值。而且我只对那些大值感兴趣，但是大量的 0 可能是奇点的原因，但需要。

4

1 回答 1

2

由于您没有提供工作示例，因此我无法轻松对此进行测试，因此举证责任在您身上。:) 无论如何，它应该是进一步修补的起点。代码中的注释。

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)
}
于 2015-12-17T10:42:03.877 回答