Как попытаться определить местонахождение ошибки, когда я использую встроенную функцию, в R? - PullRequest
Новая
       53

Как попытаться определить местонахождение ошибки, когда я использую встроенную функцию, в R?

0 голосов
/ 27 мая 2018

Я написал код в R (см. Ниже).Это работает, когда N=100.Мне нужно запустить функцию dist_statistic N=1000 раза.

Внутри этой функции неявно используется разложение Холецкого.Для разложения Холецкого матрица должна быть положительно определенной.Но элементы i -ой матрицы являются случайными числами.Я не контролирую позитивность.В результате я вижу ошибку: 

# Error in chol.default(rxx) :
# the leading minor of order 4 is not positive definite

и затем вычисления останавливаются.

Вопрос: Как отследить местоположение ошибки и продолжить вычисления с генерациейновой положительно определенной матрицы? 

library(fungible)
n <- 4
k <- 2 
p <- n 
n1 <- 100; n2 <- 100

R1 <- matrix(c(
1.00, 0.51, 0.44, 0.22,
0.51, 1.00, 0.36, 0.21,
0.44, 0.36, 1.00, 0.26,
0.22, 0.21, 0.26, 1.00), n, n)

skew_vec = c(-0.254, -0.083, 0.443, -0.017); kurt_vec = c(6.133,   4.709, 6.619,  4.276)

dist_statistic <- function(N, n, n1, n2, R1){
Q <- c()
for(i in 1:N)
{
    X1 <- monte1(seed = i+123, nvar = n, nsub = n1, cormat = R1,
             skewvec = skew_vec,
             kurtvec = kurt_vec)$data #; X1


    R2 <- corSample(R1, n = 10000)$cor.sample

    rand_vec <- rnorm(n)

    X2 <- monte1(seed = i+321, nvar = n, nsub = n2, cormat = R2,
                 skewvec = skew_vec + rand_vec,
                 kurtvec = kurt_vec + rand_vec)$data

    G1 <- adfCor(X1);    G2 <- adfCor(X2)

    G    <- ((n1 - 1)*G1 + (n2 - 1)*G2)/(n1 + n2 - 2)
    Ginv <- MASS::ginv(G)

    # vectorization operator
    delta <- row(R1) - col(R2)
    vR1 <- as.vector(t(R1[delta > 0])); vR2 <- as.vector(t(R2[delta > 0]))

    stat  <- n1*n2/(n1 + n2) * ((vR1 - vR2) %*% Ginv) %*% (vR1 - vR2)
    Q <- c(Q, stat) 
    print(i)
 } # for_i

      Results <- list(statistic = Q, iteration = i)
    return(Results)
} # function

s <- dist_statistic(N=100, n, n1, n2, R1)

1 Ответ

0 голосов
/ 27 мая 2018

Вот подход.Сначала я переписываю содержимое вашего цикла как функцию: 

my_function <- function(i) {
  X1 <- monte1(seed = i+123, nvar = n, nsub = n1, cormat = R1,
               skewvec = skew_vec,
               kurtvec = kurt_vec)$data #; X1
  R2 <- corSample(R1, n = 10000)$cor.sample
  rand_vec <- rnorm(n)
  X2 <- monte1(seed = i+321, nvar = n, nsub = n2, cormat = R2,
               skewvec = skew_vec + rand_vec,
               kurtvec = kurt_vec + rand_vec)$data
  G1 <- adfCor(X1)
  G2 <- adfCor(X2)
  G    <- ((n1 - 1)*G1 + (n2 - 1)*G2)/(n1 + n2 - 2)
  Ginv <- MASS::ginv(G)
  # vectorization operator
  delta <- row(R1) - col(R2)
  vR1 <- as.vector(t(R1[delta > 0]))
  vR2 <- as.vector(t(R2[delta > 0]))
  stat  <- n1*n2/(n1 + n2) * ((vR1 - vR2) %*% Ginv) %*% (vR1 - vR2)
  return(stat)
}

Теперь мы можем использовать эту функцию в tryCatch: 

dist_statistic <- function(N, n, n1, n2, R1){
  Q <- c()
  counter <- 1
  i <- 1
  while (counter <= N) {
    tryCatch({
      Q <- c(Q, my_function(i))
      cat(".")      
      counter <- counter + 1
    },
      error = function(e) {
        cat("*")
      },
      finally = {
        if (i %% 20 == 0) cat("\n")
        i <- i + 1
      }
  )}
  cat("\n")
  Results <- list(statistic = Q, iteration = i - 1)
  return(Results)
}

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

s <- dist_statistic(N=110, n, n1, n2, R1)
# ....................
# ....................
# ....................
# ....................
# ....................
# .*..*.......

str(s)
# List of 2
#  $ statistic: num [1:110] 5.91 2.59 5.49 5.01 1.65 ...
#  $ iteration: num 112
...