Qual é a maneira correta / padrão de verificar se a diferença é menor que a precisão da máquina?

Frequentemente, acabo em situações em que é necessário verificar se a diferença obtida está acima da precisão da máquina. Parece que para esta finalidade R tem uma variável útil: .Machine$double.eps. No entanto, quando recorro ao código-fonte R para obter orientações sobre o uso desse valor, vejo vários padrões diferentes.

Exemplos

Aqui estão alguns exemplos da statsbiblioteca:

t.test.R

if(stderr < 10 *.Machine$double.eps * abs(mx))

chisq.test.R

if(abs(sum(p)-1) > sqrt(.Machine$double.eps))

integrar.R

rel.tol < max(50*.Machine$double.eps, 0.5e-28)

lm.influence.R

e[abs(e) < 100 * .Machine$double.eps * median(abs(e))] <- 0

princomp.R

if (any(ev[neg] < - 9 * .Machine$double.eps * ev[1L]))

etc.

Questões

1. Como se pode compreender o raciocínio por trás de todos esses diferentes 10 *, 100 *, 50 *e sqrt()modificadores?
2. Existem diretrizes sobre como usar .Machine$double.epspara ajustar diferenças devido a problemas de precisão?

A precisão da máquina doubledepende de seu valor atual. .Machine$double.epsfornece a precisão quando os valores são 1. Você pode usar a função C nextAfterpara obter a precisão da máquina para outros valores.

library(Rcpp)
cppFunction("double getPrec(double x) {
  return nextafter(x, std::numeric_limits<double>::infinity()) - x;}")

(pr <- getPrec(1))
#[1] 2.220446e-16
1 + pr == 1
#[1] FALSE
1 + pr/2 == 1
#[1] TRUE
1 + (pr/2 + getPrec(pr/2)) == 1
#[1] FALSE
1 + pr/2 + pr/2 == 1
#[1] TRUE
pr/2 + pr/2 + 1 == 1
#[1] FALSE

A adição de valor aa valor bnão será alterada bquando ahouver <= metade da precisão da máquina. Verificando se a diferença é menor que a precisão da máquina <. Os modificadores podem considerar casos típicos com que frequência uma adição não mostrou uma alteração.

Em R, a precisão da máquina pode ser estimada com:

getPrecR <- function(x) {
  y <- log2(pmax(.Machine$double.xmin, abs(x)))
  ifelse(x < 0 & floor(y) == y, 2^(y-1), 2^floor(y)) * .Machine$double.eps
}
getPrecR(1)
#[1] 2.220446e-16

Cada doublevalor está representando um intervalo. Para uma adição simples, o intervalo do resultado depende da reorganização de cada somamand e também do intervalo da soma.

library(Rcpp)
cppFunction("std::vector<double> getRange(double x) {return std::vector<double>{
   (nextafter(x, -std::numeric_limits<double>::infinity()) - x)/2.
 , (nextafter(x, std::numeric_limits<double>::infinity()) - x)/2.};}")

x <- 2^54 - 2
getRange(x)
#[1] -1  1
y <- 4.1
getRange(y)
#[1] -4.440892e-16  4.440892e-16
z <- x + y
getRange(z)
#[1] -2  2
z - x - y #Should be 0
#[1] 1.9

2^54 - 2.9 + 4.1 - (2^54 + 5.9) #Should be -4.7
#[1] 0
2^54 - 2.9 == 2^54 - 2      #Gain 0.9
2^54 - 2 + 4.1 == 2^54 + 4  #Gain 1.9
2^54 + 5.9 == 2^54 + 4      #Gain 1.9

Para maior precissão Rmpfrpode ser usado.

library(Rmpfr)
mpfr("2", 1024L)^54 - 2.9 + 4.1 - (mpfr("2", 1024L)^54 + 5.9)
#[1] -4.700000000000000621724893790087662637233734130859375

Caso pudesse ser convertido em número inteiro, gmppoderia ser usado (o que está em Rmpfr).

library(gmp)
as.bigz("2")^54 * 10 - 29 + 41 - (as.bigz("2")^54 * 10 + 59)
#[1] -47

Definição de machine.eps: é o valor mais baixo  eps para o qual  1+eps não é 1

Como regra geral (assumindo uma representação de ponto flutuante com base 2):
Isso epsfaz a diferença para o intervalo 1 .. 2,
para o intervalo 2 .. 4 a precisão é 2*eps
e assim por diante.

Infelizmente, não existe uma boa regra geral aqui. É totalmente determinado pelas necessidades do seu programa.

Em R, temos all.equal como uma maneira integrada de testar a igualdade aproximada. Então você pode usar algo como (x<y) | all.equal(x,y)

i <- 0.1
 i <- i + 0.05
 i
if(isTRUE(all.equal(i, .15))) { #code was getting sloppy &went to multiple lines
    cat("i equals 0.15\n") 
} else {
    cat("i does not equal 0.15\n")
}
#i equals 0.15

O mock do Google tem vários marcadores de ponto flutuante para comparações de precisão dupla, incluindo DoubleEqe DoubleNear. Você pode usá-los em um combinador de matriz como este:

ASSERT_THAT(vec, ElementsAre(DoubleEq(0.1), DoubleEq(0.2)));

Atualizar:

As Receitas Numéricas fornecem uma derivação para demonstrar que o uso de um quociente de diferença unilateral sqrté uma boa opção de tamanho de etapa para aproximações de diferenças finitas de derivadas.

O site de artigo da Wikipedia Numerical Recipes, 3ª edição, Seção 5.7, que é as páginas 229-230 (um número limitado de visualizações de página está disponível em http://www.nrbook.com/empanel/ ).

all.equal(target, current,
           tolerance = .Machine$double.eps ^ 0.5, scale = NULL,
           ..., check.attributes = TRUE)

Essa aritmética de ponto flutuante IEEE é uma limitação bem conhecida da aritmética computacional e é discutida em vários locais:

• O FAQ no R tem uma pergunta inteira dedicada a ele: R FAQ 7.31
  • O R Inferno de Patrick Burns dedicou o primeiro "Círculo" a este problema (página 9 em diante)
  • Problema de prova de soma aritmética solicitado no meta-matemática
  • David Goldberg, "O que todo cientista da computação deve saber sobre aritmética de ponto flutuante", ACM Computing Surveys 23 , 1 (1991-03), 5-48 doi> 10.1145 / 103162.103163 ( revisão também disponível )
  • O guia de ponto flutuante - O que todo programador deve saber sobre aritmética de ponto flutuante
  • 0.30000000000000004.com compara aritmética de ponto flutuante entre linguagens de programação
• A duplicata canônica para "ponto flutuante é imprecisa" (uma meta-discussão sobre uma resposta canônica para esse problema)
• Várias questões de estouro de pilha, incluindo
  • Por que os números decimais não podem ser representados exatamente em binário?
  • Por que os números de ponto flutuante são imprecisos?
  • A matemática do ponto flutuante está quebrada?
• Truques de matemática explicados por Arthur T. Benjamin

. dplyr::near()é outra opção para testar se dois vetores de números de ponto flutuante são iguais.

A função possui um parâmetro de tolerância embutido: tol = .Machine$double.eps^0.5que pode ser ajustado. O parâmetro padrão é o mesmo que o padrão para all.equal().