Как вы можете обойти ошибку единственной матрицы, чтобы продолжить повторную выборку в цикле / репликации? - PullRequest
Новая
       23

Как вы можете обойти ошибку единственной матрицы, чтобы продолжить повторную выборку в цикле / репликации?

0 голосов
/ 29 мая 2019

Я запускаю симуляцию из двух частей на основе модели препятствий. Однако при начальной загрузке и повторной выборке я сталкиваюсь с ошибкой 

Error in solve.default(as.matrix(fit_zero$hessian)) : Lapack routine dgesv: system is exactly singular: U[1,1] = 0

После запуска traceback() меня встречают с этим: 

19: solve(as.matrix(fit_zero$hessian))
18: hurdle(formula = y ~ m + x, data = boot.data0, dist = "poisson", 
        zero.dist = "binomial") at #21
17: bootstrap(v)
16: FUN(X[[i]], ...)
15: lapply(X = X, FUN = FUN, ...)
14: sapply(integer(n), eval.parent(substitute(function(...) expr)), 
        simplify = simplify)
13: replicate(r, bootstrap(v), simplify = FALSE) at #32
12: unlist(replicate(r, bootstrap(v), simplify = FALSE)) at #32
11: matrix(unlist(replicate(r, bootstrap(v), simplify = FALSE)), 
        ncol = 8, byrow = TRUE) at #32
10: qnorm((sum(boot[, 3] > boot[, 1]) + sum(boot[, 3] == boot[, 1])/2)/r) at #3
9: bias.correct(matrix(unlist(replicate(r, bootstrap(v), simplify = FALSE)), 
       ncol = 8, byrow = TRUE), r) at #32
8: FUN(X[[i]], ...)
7: lapply(X = X, FUN = FUN, ...)
6: sapply(model_values, bootstrapper) at #35
5: FUN(newX[, i], ...)
4: apply(getParameters()[, 2:7], 1, function(parameters) {
       print(parameters)
       PROGRESS <- 0
       seed <- parameters["seed"]
       n <- parameters["n"]
       a <- parameters["a"]
       b <- parameters["b"]
       c <- parameters["c"]
       i <- parameters["i"]
       set.seed(seed)
       model_values <- replicate(iterations, model(n, a, b, c, i), 
           simplify = FALSE)
       bootstrapper <- function(v) {
           PROGRESS <<- PROGRESS + 1
           print(PROGRESS)
           bias.correct(matrix(unlist(replicate(r, bootstrap(v), 
               simplify = FALSE)), ncol = 8, byrow = TRUE), r)
       }
       boot.fit <- sapply(model_values, bootstrapper)
       boot.fit.matrix <- matrix(unlist(boot.fit), ncol = 4, byrow = TRUE)
       averaged <- apply(boot.fit.matrix, 2, mean)
       return(averaged)
   }) at #9
3: main()
2: as_mapper(.f)
1: possibly(main(), otherwise = "error")

После предыдущего, смутного вопроса о переполнении стека и большого количества исследований я понял, что распределение вполне может иметь все 0 или все 1 в этой конкретной выборке. Хотя я не могу использовать этот образец, я не хочу, чтобы цикл (или репликация в этом случае) останавливался.

Я хотел бы знать, есть ли способ создать условие, например, «если матрица сингулярна, выберите следующую итерацию».

Мои случайно сгенерированные данные поступают из этой части кода: 


gen.hurdle = function(n, a, b1, b2, c1, c2, i0, i1, i2){

  x = round(rnorm(n),3)
  e = rnorm(n)
  m = round(i0 + a*x + e, 3)

  lambda = exp(i1 + b1*m + c1*x)                       # PUT REGRESSION TERMS FOR THE CONTINUUM PART HERE; KEEP exp()
  ystar = qpois(runif(n, dpois(0, lambda), 1), lambda) # Zero-TRUNCATED POISSON DIST.; THE CONTINUUM PART

  z = i2 + b2*m  + c2*x                                # PUT REGRESSION TERMS FOR THE BINARY    PART HERE
  z_p = exp(z) / (1+exp(z))                            # p(1) = 1-p(0)
  tstar = rbinom(n, 1, z_p)                            # BINOMIAL DIST.         ; THE BINARY    PART

  y= ystar*tstar                                       # TWO-PART COUNT OUTCOME

  return(cbind(x,m,y,z,z_p,tstar))
}

# Returns the base model's powers, type 1 errors, and data
model <- function(n, a, b, c, i){
  #generate random data
  data  = data.frame(gen.hurdle(n, a, b, b, c, c, i, i, i))
  data0 = data.frame(gen.hurdle(n, a, 0, 0, c, c, i, i, i))

Я пробовал try() и trycatch() безрезультатно, потому что у меня много вложенных функций. Однако я прочитал аналогичный пост SO, в котором говорилось: «Вместо использования tryCatch вы можете просто вычислить определитель матрицы с помощью функции det. Матрица является сингулярной тогда и только тогда, когда определитель равен нулю».

Мне нужна помощь с установкой этого условия.

...