为什么chisq.test在求和之前按降序对数据进行排序

Question

为什么chisq.testR 中的函数在求和之前按降序对数据进行排序？

有问题的代码是:

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

如果由于使用浮点运算而担心数值稳定性并希望使用一些易于部署的hack,我会在求和之前按递增顺序对数据进行排序,以避免在累加器中将较小的值添加到大值(以避免尽可能多地修剪结果中的最低有效位.

我看着的源代码总和,但它没有解释为什么在传递数据下降为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.

Answer 1

编辑: " Nicolas J. Higham的数值算法的准确性和稳定性 "一书指出

"当通过递归求和对非负数求和时,增量排序是最佳排序,在具有最小的先验前向误差界限的意义上."

感谢@Lamia在评论部分分享了这本书.

本书介绍了三种求和方法,如递归,插入和成对技术.每种技术都有自己的优缺点,基于与它们相关的误差界限的大小,可以通过对浮点数求和进行系统误差分析来计算.

值得注意的是,递归技术的总和结果取决于排序策略,例如增加,减少和Psum(在书中查看 - 第82页 - 第4段.另请参阅第82页底部给出的示例中的工作原理.) .

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

浮点有效数中的基数也是53,可以从中获得

.Machine$double.digits
# [1] 53

通过将此数字设置为精度位,我们可以比较基本R和mpfr()Rmpfr库的总和操作的准确性,以用于不同的排序策略.请注意,增加顺序会产生更接近于浮点识别摘要中所见结果的结果,这证实了本书中的上述陈述.

我使用原始数据计算了卡方统计量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()函数中使用的原始数据值的符号(+ ve或-ve)无关.此函数不接受负值,因此通过使用此消息暂停函数将吐出错误"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

求和浮点数时的值的差异与双精度算法有关.请参阅下面的演示.当浮点数的精度被保持在使用200个位数mpfr()可用功能Rmpfr包,其中总和操作,而不管该矢量的顺序的给出相同的结果x1或x2.但是,当没有保持浮点精度时,会观察到不相等的值.

没有FP精度:

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 Precision保持:

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