浮点错误 fligner.test r 函数?

floating-point error fligner.test r function?

提问人:Jinglestar 提问时间:8/23/2022 最后编辑:Martin MächlerJinglestar 更新时间:8/29/2022 访问量:63

问:

我注意到,使用 r 包中的统计测试可以通过简单的转换提供不同的结果,即使情况并非如此。fligner.teststats

下面是一个示例(原始数据集的差异要大得多):

g  <- factor(rep(1:2, each=6))
x1 <- c(2,2,6,6,1,4,5,3,5,6,5,5)
x2 <- (x1-1)/5 #> cor(x1,x2) [1] 1
fligner.test(x1,g) # chi-squared = 4.2794, df = 1, p-value = 0.03858
fligner.test(x2,g) # chi-squared = 4.8148, df = 1, p-value = 0.02822

查看函数代码,我注意到中值居中可能会导致问题:

x1 <- x1 - tapply(x1,g,median)[g]
x2 <- x2 - tapply(x2,g,median)[g]
unique(abs(x1)) # 1 3 2 0
unique(abs(x2)) # 0.2 0.6 0.4 0.2 0.0 <- repeated 0.2

这是一个已知问题,应该如何解决这种不一致问题?

r 浮点 精度 float-accuracy

评论


答:

3赞 Allan Cameron 8/23/2022 #1

我认为你的分析在这里是正确的。在您的示例中,问题最终会发生,因为除非四舍五入到小数点后 15 位。您应该提交错误报告,因为这是可以避免的。(0.8 - 0.6) == 0.2FALSE

如果你在此期间感到绝望,你可以通过在中值居中阶段应用一点点四舍五入来消除浮点不等式来适应:stats:::fligner.test.default

fligner <- function (x, g, ...) 
{
  if (is.list(x)) {
    if (length(x) < 2L) 
      stop("'x' must be a list with at least 2 elements")
    DNAME <- deparse1(substitute(x))
    x <- lapply(x, function(u) u <- u[complete.cases(u)])
    k <- length(x)
    l <- lengths(x)
    if (any(l == 0)) 
      stop("all groups must contain data")
    g <- factor(rep(1:k, l))
    x <- unlist(x)
  }
  else {
    if (length(x) != length(g)) 
      stop("'x' and 'g' must have the same length")
    DNAME <- paste(deparse1(substitute(x)), "and", 
                   deparse1(substitute(g)))
    OK <- complete.cases(x, g)
    x <- x[OK]
    g <- g[OK]
    g <- factor(g)
    k <- nlevels(g)
    if (k < 2) 
      stop("all observations are in the same group")
  }
  n <- length(x)
  if (n < 2) 
    stop("not enough observations")
  x <- round(x - tapply(x, g, median)[g], 15)
  a <- qnorm((1 + rank(abs(x))/(n + 1))/2)
  a <- a - mean(a)
  v <- sum(a^2)/(n - 1)
  a <- split(a, g)
  STATISTIC <- sum(lengths(a) * vapply(a, mean, 0)^2)/v
  PARAMETER <- k - 1
  PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
  names(STATISTIC) <- "Fligner-Killeen:med chi-squared"
  names(PARAMETER) <- "df"
  METHOD <- "Fligner-Killeen test of homogeneity of variances"
  RVAL <- list(statistic = STATISTIC, parameter = PARAMETER, 
               p.value = PVAL, method = METHOD, data.name = DNAME)
  class(RVAL) <- "htest"
  return(RVAL)
}

现在,这将返回两个向量的正确结果:

fligner(x1,g)
#> 
#> Fligner-Killeen test of homogeneity of variances
#> 
#> data:  x1 and g
#> Fligner-Killeen:med chi-squared = 4.2794, df = 1, p-value = 0.03858

fligner(x2,g) 
#> 
#> Fligner-Killeen test of homogeneity of variances
#> 
#> data:  x2 and g
#> Fligner-Killeen:med chi-squared = 4.2794, df = 1, p-value = 0.03858