为什么在进行卡方检验前要将数据按降序排序再求和?

22

chisq.test函数在R中为什么要在求和之前按降序排序数据?

有关的代码如下:

STATISTIC <- sum(sort((x - E)^2/E, decreasing = TRUE))

如果我因使用浮点运算而担心数值稳定性,并希望使用一些易于部署的技巧,那么在求和之前按增加顺序对数据进行排序,以避免将微小值添加到累加器中的大值(为了尽可能减少结果中最不显着位的截断)。

我查看了sum()函数的源代码(链接),但它没有解释为什么要像传递数据时那样按降序排列。我错过了什么吗?

一个例子:

x = matrix(1.1, 10001, 1)
x[1] = 10^16   # We have a vector with 10000*1.1 and 1*10^16
c(sum(sort(x, decreasing = TRUE)), sum(sort(x, decreasing = FALSE)))

结果:

10000000000010996 10000000000011000

当我们按升序对数据进行排序时,我们得到了正确的结果。如果我们按降序对数据进行排序,则会得到一个偏差为4的结果。


3
我发现你的问题很有趣,所以我查了一些资料。看起来,在求和时按递减顺序可以在存在大量相消的情况下提高精度,即总和的绝对值<<绝对值之和。参见书籍Accuracy and stability of numerical algorithms p.82-83。但在这里,由于所有被求和的元素都是正数,因此按递增顺序求和应该更合适。 - Lamia
@Lamia;不错的链接。似乎第82页底部提到,如果数字是“浮点数太大”,那么减少可能是最好的选择。也许是因为(x - E)^2/E这样的项很大,所以做出了这个决定? - user2957945
1
@user2957945 我不认为这是正确的,因为你提到的例子中有一个非常大的浮点数M,它是{1,M,2M,-3M}的总和,这是一个严重的抵消例子,而在这里所有数字都是正数。他们还提到对于非负数的求和,递增顺序具有更好的精度(第82页和表4.1第89页)。所以我不确定为什么这里是递减顺序...也许有更具体知识的人可以发表意见? - Lamia
我的最佳猜测(也许我们永远无法弄清楚)是,这是由于过多的咖啡因燃料导致的编码深夜错误。本意可能是使用排序来避免丢失LSB,但不小心搞反了。而这个错误是如此无害,以至于直到现在,当一个好奇和敏锐的人开始尝试跟踪源代码时,才被发现。 - dww
我猜测R是一种弱类型的伪语言。 - wildplasser
1个回答

7

编辑: 尼古拉斯·J·海曼(Nicolas J. Higham)的书《数值算法的准确性和稳定性》指出:

“在通过递归求和非负数时,按递增顺序是最佳顺序,因为它具有最小的先验前向误差界。”

感谢@Lamia在评论中分享了这本书。

此书解释了三种求和方法,如递归、插入和成对技术。每种技术都有其自身的优点和缺点,基于与它们相关的误差边界的大小可以通过对浮点数求和进行系统误差分析来计算。

显然,递归技术的求和结果取决于排序策略,例如递增、递减和Psum(查看书-第82页-第4段。还可以查看第82页底部给出的示例,了解其工作原理。)。

查看R源代码中的sum()函数,该函数可以从summary.c中获取,可以得知R在其sum()函数中实现了递归方法。

此外,浮点尾数中的基本数字位数为53,这可以从下面的内容中获得:

.Machine$double.digits
# [1] 53

通过设置这个数作为精度位,我们可以比较使用基数R和 Rmpfr库中的 mpfr() 对于不同排序策略的精度加法操作的准确性。请注意,增加排序会产生结果更接近于浮点感知求和中看到的结果,这验证了本书中的上述说法。
我使用原始数据 x 计算了卡方统计量。
library('data.table')
library('Rmpfr')
x1 = matrix(c( 10^16, rep(1.1, 10000)), 
            nrow = 10001, ncol = 1)
df1 <- data.frame(x = x1)
setDT(df1)
df1[, p := rep(1/length(x), length(x))]
s_x <- df1[, sum(x)]
df1[, E := s_x * p]
df1[, chi := ((x - E)^2/E)]

precBits <- .Machine$double.digits
x_chi <- data.frame( names = c("x_asc", "x_desc", "x_fp_asc", "x_fp_desc",
                               "chi_asc", "chi_desc", "chi_fp_asc", "chi_fp_desc"))
x_chi$vals <- c( ## x
  df1[order(x), format( sum(x), digits = 22)],
  df1[order(-x), format( sum(x), digits = 22)],
  df1[order(x), format( sum(mpfr(x, precBits = precBits)), digits = 22)],
  df1[order(-x), format( sum(mpfr(x, precBits = precBits)), digits = 22)],
  ## chi
  df1[order(chi), format( sum(chi), digits = 22)],
  df1[order(-chi), format( sum(chi), digits = 22)],
  df1[order(chi), format( sum(mpfr(chi, precBits = precBits)), digits = 22)],
  df1[order(-chi), format( sum(mpfr(chi, precBits = precBits)), digits = 22)])

x_chi
#         names                    vals
# 1       x_asc       10000000000011000
# 2      x_desc       10000000000010996
# 3    x_fp_asc 10000000000011000.00000
# 4   x_fp_desc 10000000000020000.00000
# 5     chi_asc    99999999999890014218
# 6    chi_desc    99999999999890030592
# 7  chi_fp_asc 99999999999890014208.00
# 8 chi_fp_desc 99999999999833554944.00

查看edit(chisq.test)函数的源代码可知,其中没有涉及排序操作。

此外,正如评论部分所指出的那样,它与chisq.test()函数中原始数据的值的符号(+或-)无关。此函数不接受负值,因此会通过停止函数并显示以下错误信息来抛出错误:"all entries of 'x' must be nonnegative and finite"

set.seed(2L)
chisq.test(c(rnorm(10, 0, 1)))
# Error in chisq.test(c(rnorm(10, 0, 1))) : 
#   all entries of 'x' must be nonnegative and finite

涉及到浮点数求和时差值是由于双精度算术所造成的。请参见下面的演示。使用Rmpfr包中可用的mpfr()函数将浮点数的精度保持在200位数字,向量x1x2的顺序不同,求和操作会给出相同的结果。然而,当没有保持浮点精度时,会观察到不相等的值。

无浮点精度:

x1 = matrix(c( 10^16, rep(1.1, 10000)), 
            nrow = 10001, ncol = 1)
## reverse
x2 = matrix(c( rep(1.1, 10000), 10^16 ), 
            nrow = 10001, ncol = 1)

c( format(sum(x1), digits = 22), 
   format(sum(x2), digits = 22))
# [1] "10000000000010996" "10000000000011000"

FP精度保持不变:

library('Rmpfr')
##
prec <- 200
x1 = matrix(c( mpfr( 10^16, precBits = prec),
              rep( mpfr(1.1, precBits = prec), 10000)), 
           nrow = 10001, ncol = 1)

## reverse
x2 = matrix(c( rep(mpfr(1.1, precBits = prec), 10000), 
              mpfr( 10^16, precBits = prec) ), 
           nrow = 10001, ncol = 1)
c( sum(x1), sum(x2))
# 2 'mpfr' numbers of precision  200   bits 
# [1] 10000000000011000.000000000000888178419700125232338905334472656
# [2] 10000000000011000.000000000000888178419700125232338905334472656

在基数 R 中,最小的正浮点数可以从以下代码中获得,任何小于这个数的数字都将在基数 R 中截断,这会导致求和操作的结果不同。

.Machine$double.eps
# [1] 2.220446e-16

chisq.test()函数中双精度算法感知和不感知函数的比较。

chisq.test()函数中提取相关部分,并制作一个新的函数chisq.test2()。在内部,您将看到使用mpfr()函数应用250位双精度感知前后进行比较的选项,以计算卡方统计量。可以看到,对于浮点数感知函数,结果相同,但对于原始数据则不同。

# modified chi square function:
chisq.test2 <- function (x, precBits) 
{
  if (is.matrix(x)) {
    if (min(dim(x)) == 1L) 
      x <- as.vector(x)
  }

  #before fp precision
  p = rep(1/length(x), length(x))
  n <- sum(x)
  E <- n * p

  # after fp precision
  x1 <- mpfr(x, precBits = precBits)
  p1 = rep(1/length(x1), length(x1))
  n1 <- sum(x1)
  E1 <- n1 * p1

  # chisquare statistic
  STATISTIC <- c(format(sum((x - E)^2/E), digits=22),           # before decreasing
                 format(sum(sort((x - E)^2/E, decreasing = FALSE)), digits=22), # before increasing
                 sum((x1 - E1)^2/E1),                           # after decreasing 
                 sum(sort((x1 - E1)^2/E1, decreasing = FALSE))) # after increasing

  return(STATISTIC)
}

# data
x1 = matrix(c( 10^16, rep(1.1, 10000)), 
            nrow = 10001, ncol = 1)

chisq.test2(x = x1, precBits=250)

输出:

# [[1]]  # before fp decreasing
# [1] "99999999999890030592"
# 
# [[2]]  # before fp increasing
# [1] "99999999999890014218"
# 
# [[3]]  # after fp decreasing 
# 'mpfr1' 99999999999889972569.502489584522352514811399898444554440067408531548230046685
# 
# [[4]]  # after fp increasing
# 'mpfr1' 99999999999889972569.502489584522352514811399898444554440067408531548230251906

1
你知道为什么第57行代码 STATISTIC <- sum(sort((x - E)^2/E, decreasing = TRUE)) 要进行递减排序吗?谢谢。 - user2957945
2
我有什么遗漏,还是这并没有真正回答问题? - dww
1
为什么要对这个答案进行负评?提问者使用了一个列的矩阵数据。请查看chisq.test()函数的源代码。它从未使用源代码中第57行,@user2957945在上面的评论中提到的那一行。相反,问题中给出的数据使用了第94行。 - Sathish
另外需要注意的是,按照降序排序数据与按升序排序数据相比,会导致更不精确的数值。在没有浮点精度意识的情况下,参考值为10000000000011000并按升序排列,接近于从FP精度感知数据之和获得的值。 - Sathish
1
@user2957945 我不知道为什么这个函数的开发者采用了不太精确的方法。我也不知道该如何猜测... 如果你看到代码块在第57行(pearson chisq)和第94行(给定概率的chisq),两种方法的区别在于计算pvalue的方式,即使用统计量来判断统计量在chisq分布中是否异常出现。 - Sathish

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接