Quelle est la manière correcte / standard de vérifier si la différence est inférieure à la précision de la machine?

Je me retrouve souvent dans des situations où il est nécessaire de vérifier si la différence obtenue est supérieure à la précision de la machine. On dirait à cette fin R a une variable à portée de main: .Machine$double.eps. Cependant, lorsque je me tourne vers le code source R pour obtenir des instructions sur l'utilisation de cette valeur, je vois plusieurs modèles différents.

Exemples

Voici quelques exemples tirés de la statsbibliothèque:

t.test.R

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

chisq.test.R

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

integrer.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.

Des questions

1. Comment peut - on comprendre le raisonnement derrière tous ces différents 10 *, 100 *, 50 *et sqrt()modificateurs?
2. Existe-t-il des directives sur l'utilisation .Machine$double.epspour ajuster les différences en raison de problèmes de précision?

La précision de la machine doubledépend de sa valeur actuelle. .Machine$double.epsdonne la précision lorsque les valeurs sont 1. Vous pouvez utiliser la fonction C nextAfterpour obtenir la précision de la machine pour d'autres valeurs.

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

La valeur ajoutée aà la valeur bne changera pas bquand aest la <= moitié de la précision de la machine. Vérifier si la différence est plus petite que la précision de la machine <. Les modificateurs peuvent prendre en compte les cas typiques de la fréquence à laquelle un ajout n'a pas montré de changement.

Dans R, la précision de la machine peut être estimée avec:

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

Chaque doublevaleur représente une plage. Pour un ajout simple, la plage du résultat dépend de la réanimation de chaque somme et également de la plage de leur somme.

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

Pour une plus grande précision, on Rmpfrpourrait utiliser.

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

Au cas où il pourrait être converti en entier gmppourrait être utilisé (ce qui est dans Rmpfr).

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

Définition d'une machine.eps: c'est la valeur la plus basse  eps pour laquelle  1+eps n'est pas 1

En règle générale (en supposant une représentation en virgule flottante avec la base 2):
Cela epsfait la différence pour la plage 1 .. 2,
pour la plage 2 .. 4 la précision est 2*eps
et ainsi de suite.

Malheureusement, il n'y a pas de bonne règle empirique ici. Il est entièrement déterminé par les besoins de votre programme.

Dans R, nous avons all.equal comme méthode intégrée pour tester l'égalité approximative. Vous pourriez donc utiliser quelque chose comme (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

Google Mock dispose d'un certain nombre de comparateurs à virgule flottante pour les comparaisons à double précision, y compris DoubleEqet DoubleNear. Vous pouvez les utiliser dans un matcher de tableau comme celui-ci:

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

Mise à jour:

Les recettes numériques fournissent une dérivation pour démontrer que l'utilisation d'un quotient de différence unilatérale sqrtest un bon choix de taille de pas pour les approximations de différences finies de dérivés.

Le site d'articles Wikipédia Recettes numériques, 3e édition, section 5.7, qui se trouve aux pages 229-230 (un nombre limité de pages vues est disponible sur http://www.nrbook.com/empanel/ ).

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

Cette arithmétique à virgule flottante IEEE est une limitation bien connue de l'arithmétique informatique et est discutée à plusieurs endroits:

• La FAQ de R a toute une question qui lui est dédiée: R FAQ 7.31
  • Le R Inferno de Patrick Burns a dédié le premier "Circle" à ce problème (page 9 et suivantes)
  • Problème de preuve de somme arithmétique posé dans Math Meta
  • David Goldberg, «Ce que tout informaticien devrait savoir sur l'arithmétique à virgule flottante», ACM Computing Surveys 23 , 1 (1991-03), 5-48 doi> 10.1145 / 103162.103163 ( révision également disponible )
  • Le guide à virgule flottante - Ce que chaque programmeur devrait savoir sur l'arithmétique à virgule flottante
  • 0.30000000000000004.com compare l'arithmétique à virgule flottante à travers les langages de programmation
• Duplication canonique pour "virgule flottante est inexacte" (une méta-discussion sur une réponse canonique pour ce problème)
• Plusieurs questions sur le débordement de pile, y compris
  • Pourquoi les nombres décimaux ne peuvent-ils pas être représentés exactement en binaire?
  • Pourquoi les nombres à virgule flottante sont-ils inexacts?
  • Les mathématiques en virgule flottante sont-elles cassées?
• Astuces mathématiques expliquées par Arthur T. Benjamin

. dplyr::near()est une autre option pour tester si deux vecteurs de nombres à virgule flottante sont égaux.

La fonction a un paramètre de tolérance intégré: tol = .Machine$double.eps^0.5qui peut être ajusté. Le paramètre par défaut est le même que celui par défaut pour all.equal().