如何修复R语言中while循环偶发无限运行的问题？

第一个while循环（Phase I异常点移除）的核心问题 #

原代码存在两个致命逻辑错误：

• 删除异常点后未重新计算T2All，导致Compare始终基于初始观测数据，无法正确识别新数据中的异常，进而引发无限循环
• 循环条件isTRUE(Compare > UCL)完全错误：Compare是异常点的索引向量，和控制限数值UCL没有可比性；当无异常点时Compare为空向量，isTRUE()无法正确判断终止条件

修复后的Phase I循环代码 #

# 修正：循环检测并移除异常点，直到无异常点为止
while(length(Compare) > 0) {
    # 删除异常点
    Obs <- Obs[-Compare,]
    
    # 重新计算均值、协方差和T2统计量
    VecT2 <- apply(Obs, 2, mean)
    VecT2 <- round(VecT2, 3)
    
    ST2 <- cov(Obs)
    ST2 <- round(ST2, 3)
    
    T2All <- apply(Obs, 1, function(x) t(x - VecT2) %*% solve(ST2) %*% (x - VecT2))
    
    # 重新计算控制限和异常点索引
    Alpha <- 0.005
    M <- nrow(Obs)
    p <- ncol(Obs)
    
    UCL <- ((p * (M-1) * (M + 1))) / ((M - p) * M) * qf((1-Alpha), p, (M-p))
    UCL <- round(UCL, 3)
    
    Compare <- which(T2All > UCL)
}

第二个while循环（ARL0异常处理）的潜在问题 #

原代码在极端情况下可能陷入无限循环：若持续生成的Obs20_1中存在异常，循环会一直运行。此外，当Result为NA（前20个观测均未超控制限）时，原条件判断逻辑虽能终止循环，但需添加循环次数限制避免极端情况。

修复后的ARL0循环代码 #

# 修正：添加最大循环次数限制，防止无限运行
max_iter <- 1000  # 可根据需求调整
iter_count <- 0

while((!is.na(Result) && Result < 20) && iter_count < max_iter) {
    iter_count <- iter_count + 1
    
    Obs20_1 <- mvrnorm(20, mu = mu, Sigma = Sigma)
    Obs40 <- rbind(Obs20_1, Obs20_2)
    Obs40 <- as.matrix(Obs40)
    
    T2 <- apply(Obs40, 1, function(x) t(x - mu) %*% solve(Sigma) %*% (x - mu))
    Result <- which(T2 > UCL)[1]
}

# 达到最大循环次数时给出提示
if(iter_count == max_iter) {
    warning("达到最大循环次数，可能存在异常情况")
}

其他优化建议 #

• 移除冗余代码：mu <- t(mu); mu <- t(mu) 可直接替换为mu <- matrix(c(0,0,0), ncol=1)
• 用向量化操作替代for循环计算T2统计量，提升运行效率
• 移除未使用的包：readxl、dplyr、ggplot2在代码中未发挥作用，可删除

完整修复后的代码 #

rm(list=ls())

library(MASS)

# Mean Vector, Covariance Matrix Construction
mu <- matrix(c(0,0,0), ncol=1)
mu2 <- matrix(c(1,2,1), ncol=1)
Sigma <- matrix(c(1, 0.9, 0.9, 0.9, 1, 0.9, 0.9, 0.9, 1), 3)

getResult <- function() {
    # Construct 50 Random Variables for Phase I
    Obs <- mvrnorm(50, mu = mu, Sigma = Sigma)
    
    VecT2 <- apply(Obs, 2, mean)
    VecT2 <- round(VecT2, 3)
    
    ST2 <- cov(Obs)
    ST2 <- round(ST2, 3)
    
    Obs <- as.matrix(Obs)    
    T2All <- apply(Obs, 1, function(x) t(x - VecT2) %*% solve(ST2) %*% (x - VecT2))
    
    # Construct Control Limit
    Alpha <- 0.005
    M <- nrow(Obs)
    p <- ncol(Obs)
    UCL <- ((p * (M-1) * (M + 1))) / ((M - p) * M) * qf((1-Alpha), p, (M-p))
    UCL <- round(UCL, 3)
    
    Compare <- which(T2All > UCL)
    
    # Repeat removing out-of-control points in Phase I
    while(length(Compare) > 0) {
        Obs <- Obs[-Compare,]
        
        # Recalculate statistics with updated Obs
        VecT2 <- apply(Obs, 2, mean)
        VecT2 <- round(VecT2, 3)
        
        ST2 <- cov(Obs)
        ST2 <- round(ST2, 3)
        
        T2All <- apply(Obs, 1, function(x) t(x - VecT2) %*% solve(ST2) %*% (x - VecT2))
        
        # Recalculate UCL and Compare
        M <- nrow(Obs)
        UCL <- ((p * (M-1) * (M + 1))) / ((M - p) * M) * qf((1-Alpha), p, (M-p))
        UCL <- round(UCL, 3)
        
        Compare <- which(T2All > UCL)
    }
    
    # Prepare Observations for Phase II
    Obs20_1 <- mvrnorm(20, mu = mu, Sigma = Sigma)
    Obs20_2 <- mvrnorm(20, mu = mu2, Sigma = Sigma)
    Obs40 <- rbind(Obs20_1, Obs20_2)
    Obs40 <- as.matrix(Obs40)
    
    T2 <- apply(Obs40, 1, function(x) t(x - mu) %*% solve(Sigma) %*% (x - mu))
    Result <- which(T2 > UCL)[1]
    
    # Repeat until no out-of-control in first 20 observations
    max_iter <- 1000
    iter_count <- 0
    while((!is.na(Result) && Result < 20) && iter_count < max_iter) {
        iter_count <- iter_count + 1
        
        Obs20_1 <- mvrnorm(20, mu = mu, Sigma = Sigma)
        Obs40 <- rbind(Obs20_1, Obs20_2)
        Obs40 <- as.matrix(Obs40)
        
        T2 <- apply(Obs40, 1, function(x) t(x - mu) %*% solve(Sigma) %*% (x - mu))
        Result <- which(T2 > UCL)[1]
    }
    
    if(iter_count == max_iter) {
        warning("Reached maximum iterations for ARL0 check")
    }
    
    Result
}

# Run simulation
Final <- replicate(n = 200, expr = getResult())
Final <- Final - 20
print(Final)
cat("Mean of Final:", mean(Final, na.rm = TRUE), "\n")

内容的提问来源于stack exchange，提问作者Joseph Kim