编辑:
尼古拉斯·J·海曼(Nicolas J. Higham)的书《数值算法的准确性和稳定性》指出:
“在通过递归求和非负数时,按递增顺序是最佳顺序,因为它具有最小的先验前向误差界。”
感谢@Lamia在评论中分享了这本书。
此书解释了三种求和方法,如递归、插入和成对技术。每种技术都有其自身的优点和缺点,基于与它们相关的误差边界的大小可以通过对浮点数求和进行系统误差分析来计算。
显然,递归技术的求和结果取决于排序策略,例如递增、递减和Psum(查看书-第82页-第4段。还可以查看第82页底部给出的示例,了解其工作原理。)。
查看R源代码中的sum()
函数,该函数可以从summary.c中获取,可以得知R在其sum()
函数中实现了递归方法。
此外,浮点尾数中的基本数字位数为53,这可以从下面的内容中获得:
.Machine$double.digits
通过设置这个数作为精度位,我们可以比较使用基数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(
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)],
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
查看edit(chisq.test)
函数的源代码可知,其中没有涉及排序操作。
此外,正如评论部分所指出的那样,它与chisq.test()
函数中原始数据的值的符号(+或-)无关。此函数不接受负值,因此会通过停止函数并显示以下错误信息来抛出错误:"all entries of 'x' must be nonnegative and finite"
。
set.seed(2L)
chisq.test(c(rnorm(10, 0, 1)))
涉及到浮点数求和时差值是由于双精度算术所造成的。请参见下面的演示。使用
Rmpfr
包中可用的
mpfr()
函数将浮点数的精度保持在200位数字,向量
x1
或
x2
的顺序不同,求和操作会给出相同的结果。然而,当没有保持浮点精度时,会观察到不相等的值。
无浮点精度:
x1 = matrix(c( 10^16, rep(1.1, 10000)),
nrow = 10001, ncol = 1)
x2 = matrix(c( rep(1.1, 10000), 10^16 ),
nrow = 10001, ncol = 1)
c( format(sum(x1), digits = 22),
format(sum(x2), digits = 22))
FP精度保持不变:
library('Rmpfr')
prec <- 200
x1 = matrix(c( mpfr( 10^16, precBits = prec),
rep( mpfr(1.1, precBits = prec), 10000)),
nrow = 10001, ncol = 1)
x2 = matrix(c( rep(mpfr(1.1, precBits = prec), 10000),
mpfr( 10^16, precBits = prec) ),
nrow = 10001, ncol = 1)
c( sum(x1), sum(x2))
在基数 R 中,最小的正浮点数可以从以下代码中获得,任何小于这个数的数字都将在基数 R 中截断,这会导致求和操作的结果不同。
.Machine$double.eps
chisq.test()
函数中双精度算法感知和不感知函数的比较。
从chisq.test()
函数中提取相关部分,并制作一个新的函数chisq.test2()
。在内部,您将看到使用mpfr()
函数应用250位双精度感知前后进行比较的选项,以计算卡方统计量。可以看到,对于浮点数感知函数,结果相同,但对于原始数据则不同。
chisq.test2 <- function (x, precBits)
{
if (is.matrix(x)) {
if (min(dim(x)) == 1L)
x <- as.vector(x)
}
p = rep(1/length(x), length(x))
n <- sum(x)
E <- n * p
x1 <- mpfr(x, precBits = precBits)
p1 = rep(1/length(x1), length(x1))
n1 <- sum(x1)
E1 <- n1 * p1
STATISTIC <- c(format(sum((x - E)^2/E), digits=22),
format(sum(sort((x - E)^2/E, decreasing = FALSE)), digits=22),
sum((x1 - E1)^2/E1),
sum(sort((x1 - E1)^2/E1, decreasing = FALSE)))
return(STATISTIC)
}
x1 = matrix(c( 10^16, rep(1.1, 10000)),
nrow = 10001, ncol = 1)
chisq.test2(x = x1, precBits=250)
输出: