R矩阵求和浮点求和边距误差

Question

我有这个数据：

   SEX LS  GPA
1    M DO 2.56
2    M DO 2.70
3    M DO 3.04
4    M DO 1.42
5    M FS 3.02
6    M FS 2.67
7    M FS 2.35
8    M FS 1.80
9    M CO 3.41
10   M CO 2.97
11   M CO 2.23
12   M CO 2.31
13   M AP 2.28
14   M AP 3.73
15   M AP 2.05
16   M AP 2.61
17   M HO 2.57
18   M HO 2.81
19   M HO 3.01
20   M HO 3.40
21   F DO 3.33
22   F DO 1.80
23   F DO 2.50
24   F DO 3.04
25   F FS 3.87
26   F FS 3.00
27   F FS 3.25
28   F FS 2.76
29   F CO 3.14
30   F CO 4.00
31   F CO 2.66
32   F CO 2.91
33   F AP 3.69
34   F AP 2.55
35   F AP 3.21
36   F AP 2.86
37   F HO 3.09
38   F HO 1.99
39   F HO 3.46
40   F HO 3.61

我是运行数据方差分析，为了找到参数，我正在执行以下操作：

计算平均值：

mu_grand <- mean(data$GPA)
mu_alpha <- aggregate(data$GPA,list(data$SEX),mean)$x
mu_beta <- aggregate(data$GPA,list(data$LS),mean)$x
mu_alpha_beta <- matrix(aggregate(data$GPA,list(data$SEX,data$LS),mean)$x,nrow=2)

计算加性效应

alpha <- mu_alpha - mu_grand
beta <- mu_beta - mu_grand

计算交互效果

a1 <- rep(alpha[1],length(beta))
a2 <- rep(alpha[2],length(beta))

ab1 <- a1 + beta
ab2 <- a2 + beta

ab <- rbind(ab1,ab2)

alpha_beta <- mu_alpha_beta - (mu_grand + ab)

这里是alpha_beta:

ab1  0.0105  0.02925 -0.07575  0.1855 -0.1495
ab2 -0.0105 -0.02925  0.07575 -0.1855  0.1495

如果我想求 alpha_beta 的总和，我会得到

sum(alpha_beta)
[1] -4.440892e-16

虽然这绝对在浮点精度范围内，但我想知道如果您只是手动将所有数字相加得到 0，为什么会出现任何错误？谢谢！

Answer 1

你看那个数字，真是少之又少。 R 只计算到大约 15 位数字，所以这只是舍入误差。使用包 gmp 我们可以看到物理定律并没有被打破。

加载需要的包

require("gmp") #Required package

读入数据

data<- read.table(text="SEX LS  GPA
1    M DO 2.56
                  2    M DO 2.70
                  3    M DO 3.04
                  4    M DO 1.42
                  5    M FS 3.02
                  6    M FS 2.67
                  7    M FS 2.35
                  8    M FS 1.80
                  9    M CO 3.41
                  10   M CO 2.97
                  11   M CO 2.23
                  12   M CO 2.31
                  13   M AP 2.28
                  14   M AP 3.73
                  15   M AP 2.05
                  16   M AP 2.61
                  17   M HO 2.57
                  18   M HO 2.81
                  19   M HO 3.01
                  20   M HO 3.40
                  21   F DO 3.33
                  22   F DO 1.80
                  23   F DO 2.50
                  24   F DO 3.04
                  25   F FS 3.87
                  26   F FS 3.00
                  27   F FS 3.25
                  28   F FS 2.76
                  29   F CO 3.14
                  30   F CO 4.00
                  31   F CO 2.66
                  32   F CO 2.91
                  33   F AP 3.69
                  34   F AP 2.55
                  35   F AP 3.21
                  36   F AP 2.86
                  37   F HO 3.09
                  38   F HO 1.99
                  39   F HO 3.46
                  40   F HO 3.61", header=T,row.names=1)

代码

data$GPA<-as.bigq(data$GPA) #Converts to class `bigq`

#Special function because sapply was converting data to raw format
meanL<-function(L){ 
  x<-rep(0,length(L))
  x<-as.bigq(x)

  for(i in 1:length(L)){
    x[i]<-mean(L[[i]])
  }
  return(x)
}

#This section is changed to handle the new class `bigq`
mu_grand <- mean(data$GPA)
mu_alpha <- meanL(split(data$GPA,data$SEX))
mu_beta <- meanL(split(data$GPA,data$LS))
mu_alpha_beta <- matrix(meanL(split(data$GPA,as.factor(paste(data$SEX,data$LS,sep="-")))),nrow=2)



alpha <- mu_alpha - mu_grand
beta <- mu_beta - mu_grand

a1 <- rep(alpha[1],length(beta))
a2 <- rep(alpha[2],length(beta))

ab1 <- a1 + beta
ab2 <- a2 + beta

ab<-matrix(c(ab1,ab2),ncol=5,byrow=T) #The same as rbind but needed to work with objects of class 'bigq'.
alpha_beta <- mu_alpha_beta - (mu_grand + ab)
sum(alpha_beta)

Big Rational ('bigq') :
[1] 0

这是一个很好的思想实验，但它在 **s 中工作很痛苦 gmp 因为你可以看到像 aggregate 和 sapply 这样的函数不玩很好。这个故事的寓意是当你看到非常小的数字，即 1e-16 而你期待 0 你可以假设它们是 0.

R矩阵求和浮点求和边距误差

R Matrix Sum Floating Point Sum Margin Error

r

anova