22

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

有问题的代码是：

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

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

我查看了sum的源代码，但它没有解释为什么要按降序将数据传递给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 的结果。

4

1 回答 1

7

编辑：Nicolas J. Higham 所著的“数值算法的准确性和稳定性”一书指出

“当通过递归求和对非负数求和时，从具有最小先验前向误差界限的意义上说，递增排序是最好的排序。”

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

这本书解释了三种求和方法，例如递归、插入和成对技术。根据与它们相关的误差范围的大小，每种技术都有自己的优点和缺点，可以通过对浮点数求和进行系统误差分析来计算。

值得注意的是，递归技术的求和结果取决于增加、减少和 Psum 等排序策略（请参阅本书 - 第 82 页 - 第 4 段。另请参阅第 82 页底部给出的示例中它是如何工作的。） .

查看sum()可以从summary.csum()获得的函数的 R 源代码，可以发现 R 在其函数中实现了递归方法。

浮点有效数中的基数也是 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

对浮点数求和时的值差异与双精度算术有关。请参阅下面的演示。mpfr()当使用包中提供的函数将浮点数的精度保持在 200 位时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
于 2017-12-26T12:11:29.727 回答