逐步截尾均匀顺序统计量CDF大数运算异常原因及解决方案

我正在求解逐步截尾均匀顺序统计量的累积分布函数（CDF），编写的R代码如下：

#The CDF for rth order progressively censored uniform order statistics
n <- 30 #Total number of experimental units
m <- 15 #Desired numbers of failure
R <- c(rep(0, m - 1), n - m) #Progressive censoring scheme
order <- m #Order of the censored order statistics (Here, maximum order)

gam <- NA
Cr <- NA
for(i in 1 : order)
{
  gam[i] <- m - i + 1 + sum(R[i : m])
}
for(i in 1 : order)
{
  Cr[i] <- prod(gam[1 : i])
}
air <- array(dim = c(order, order))

for(i in 1 : order)
{
  for (j in 1 : order) {
    
    if(i != j)
    {
      air[i, j] <- 1/(gam[j] - gam[i])
    }
  }
}
A <- NA
for(i in 1 : order)
{
  A[i] = prod(na.omit(air[i,]))
}
#CDF of progressively censored uniform order statistics
progU_CDF <- function(u)
{
  CDF = NA
  for(i in 1 : length(u))
  {
    CDF[i] <- 1 - (Cr[order] * sum((A/gam) * ((1 - u[i])^(gam))))
  }
  return(CDF)
}

理论上，progU_CDF(0)应返回0，progU_CDF(1)应返回1。当n=30、m=15时，progU_CDF(0)的结果非常接近0，CDF曲线呈单调非减的理想形态；但当参数改为n=50、m=25时，Cr[order] * sum(A/gam)的结果远偏离1，CDF曲线形态异常。此外，Cr[order] * sum(A/gam)与sum(Cr[order] * A/gam)的计算结果不同，这与Cr[order]为常量的直觉相悖。

问题原因分析 #

1. 大数与极小值相乘的精度丢失
   当m=25、n=50时，Cr[order]是gam序列前25项的乘积：50*49*...*26，这是一个远超浮点数有效精度范围的极大数；而A[i]是1/(gam[j]-gam[i])（j≠i）的乘积，由于gam[j]-gam[i]是整数，A[i]是极小的小数（甚至接近机器精度下限）。极大数与极小值直接相乘时，浮点数的有效位数有限，会导致大量精度丢失，计算结果严重偏离理论值。

2. 浮点数运算顺序的影响
   Cr[order] * sum(A/gam)与sum(Cr[order] * A/gam)结果不同，本质是浮点数舍入误差的累积差异：先求和再乘大数，会因求和结果的精度不足导致最终误差；先乘大数再求和，会因每个乘积项的精度丢失（大数乘小数）导致误差累积，两种顺序的误差路径不同，结果自然不一致。

解决方案：基于对数的高精度计算 #

核心思路是将乘法/除法转换为对数的加法/减法，避免直接处理大数和极小值，利用R的lgamma函数（对数伽马函数，可高效计算阶乘的对数）来保持数值精度。

修改后的代码如下：

# 逐步截尾均匀顺序统计量CDF（高精度版本）
progU_CDF_high_precision <- function(n, m, u) {
  # 直接推导gam序列：n, n-1, ..., n-m+1
  gam <- n - 0:(m-1)
  
  # 计算每个项的对数绝对值和符号
  log_abs_terms <- numeric(m)
  sign_terms <- numeric(m)
  
  for (i in 1:m) {
    # Cr[m] = n!/(n-m)!，转换为对数计算
    log_Cr <- lgamma(n + 1) - lgamma(n - m + 1)
    # sum(log|i-j|) = log((i-1)!) + log((m-i)!)，用lgamma计算
    log_prod_abs_diff <- lgamma(i) + lgamma(m - i + 1)
    
    log_abs_terms[i] <- log_Cr - log(gam[i]) - log_prod_abs_diff
    # 符号部分：(-1)^(m-i)
    sign_terms[i] <- (-1)^(m - i)
  }
  
  # 计算每个项的实际值
  terms <- sign_terms * exp(log_abs_terms)
  
  # 计算CDF，处理边界值避免无效运算
  cdf <- numeric(length(u))
  for (idx in 1:length(u)) {
    u_val <- u[idx]
    if (u_val >= 1) {
      cdf[idx] <- 1
    } else if (u_val <= 0) {
      cdf[idx] <- 0
    } else {
      # 用对数计算(1-u)^gam，避免指数溢出
      log_1_minus_u_pow_gam <- gam * log(1 - u_val)
      term_contribution <- terms * exp(log_1_minus_u_pow_gam)
      cdf[idx] <- 1 - sum(term_contribution)
    }
  }
  
  return(cdf)
}

关键改进点：

• 直接推导gam序列：根据截尾方案，gam[i] = n - i + 1，无需循环计算，简化代码。
• 对数转换计算：用lgamma计算阶乘/排列数的对数，将所有乘法/除法转换为对数的加减，彻底避免大数和极小值的直接运算。
• 符号单独处理：将项的符号与绝对值分离计算，避免对数运算处理符号的问题。
• 边界值直接返回：对u<=0和u>=1直接返回理论值，避免无效计算。

验证效果 #

调用progU_CDF_high_precision(50,25,0)会返回0，progU_CDF_high_precision(50,25,1)返回1；对于任意n和m，Cr[order] * sum(A/gam)的理论值1会被精确计算，CDF曲线保持单调非减的理想形态。

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