在 R 中使用内置函数时如何尝试捕获错误位置?
How to try-and-catch error location when I use the buil-in function, in R?
我已经在R
中编写了代码(见下文)。它在 N=100
时有效。
我需要 运行 dist_statistic 函数 N=1000
次。
在此函数中,隐式使用了 Cholesky 分解。对于 Cholesky 分解,矩阵必须是正定的。但是第 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)
这是一种方法。我首先将你的循环内容重写为一个函数:
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 纯粹是为了装饰目的,表示错误。因此
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
我已经在R
中编写了代码(见下文)。它在 N=100
时有效。
我需要 运行 dist_statistic 函数 N=1000
次。
在此函数中,隐式使用了 Cholesky 分解。对于 Cholesky 分解,矩阵必须是正定的。但是第 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)
这是一种方法。我首先将你的循环内容重写为一个函数:
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 纯粹是为了装饰目的,表示错误。因此
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