Jaki jest prawidłowy / standardowy sposób sprawdzenia, czy różnica jest mniejsza niż precyzja maszyny?

Często kończę w sytuacjach, w których konieczne jest sprawdzenie, czy uzyskana różnica jest wyższa niż precyzja maszyny. Wygląda na to, w tym celu R ma zmienną poręczny: .Machine$double.eps. Jednak po przejściu do kodu źródłowego R w celu uzyskania wskazówek dotyczących korzystania z tej wartości widzę wiele różnych wzorców.

Przykłady

Oto kilka przykładów z statsbiblioteki:

t.test.R

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

chisq.test.R

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

integrate.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]))

itp.

pytania

1. Jak można zrozumieć uzasadnienie wszystkich tych różnych 10 *, 100 *, 50 *oraz sqrt()modyfikatory?
2. Czy istnieją wytyczne dotyczące używania .Machine$double.epsdo korygowania różnic z powodu problemów z precyzją?

Precyzja maszyny doublezależy od jej aktualnej wartości. .Machine$double.epsdaje precyzję, gdy wartości wynoszą 1. Możesz użyć funkcji C, nextAfteraby uzyskać dokładność maszyny dla innych wartości.

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

Dodanie wartości ado wartości bnie zmieni się, bgdy abędzie to <= połowa jej precyzji maszyny. Dokonuje się sprawdzenia, czy różnica jest mniejsza niż precyzja maszyny <. Modyfikatory mogą brać pod uwagę typowe przypadki, jak często dodatek nie pokazywał zmiany.

W R precyzję maszyny można oszacować za pomocą:

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

Każda doublewartość reprezentuje zakres. Dla prostego dodania, zakres wyniku zależy od zmiany każdego summanda, a także zakresu ich sumy.

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

Dla wyższej precyzji Rmpfrmożna zastosować.

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

W przypadku, gdy można go przekonwertować na liczbę całkowitą, gmpmożna użyć (co jest w Rmpfr).

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

Definicja machine.eps: jest to najniższa wartość,  eps dla której  1+eps nie jest 1

Zasadniczo (przy założeniu reprezentacji zmiennoprzecinkowej z bazą 2):
To epsrobi różnicę dla zakresu 1 .. 2,
dla zakresu 2 .. 4 precyzja jest 2*eps
i tak dalej.

Niestety nie ma tutaj dobrej zasady. Jest to całkowicie zależne od potrzeb Twojego programu.

W R mamy all.equal jako wbudowany sposób testowania przybliżonej równości. Więc możesz użyć czegoś takiego (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 ma wiele dopasowań zmiennoprzecinkowych do porównań podwójnej precyzji, w tym DoubleEqi DoubleNear. Możesz użyć ich w dopasowaniu tablic, takim jak to:

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

Aktualizacja:

Przepisy numeryczne dostarczają pochodnych, aby wykazać, że zastosowanie jednostronnego ilorazu różnicy sqrtjest dobrym wyborem wielkości kroku dla przybliżeń różnic skończonych pochodnych.

Witryna z artykułami Wikipedii Przepisy numeryczne, wydanie trzecie, sekcja 5.7, czyli strony 229–230 (ograniczona liczba odsłon jest dostępna pod adresem http://www.nrbook.com/empanel/ ).

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

Te arytmetyki zmiennoprzecinkowe IEEE są dobrze znanym ograniczeniem arytmetyki komputerowej i są omawiane w kilku miejscach:

• FAQ w R ma na to całe pytanie: R FAQ 7.31
  • R Inferno Patricka Burnsa poświęcił temu problemowi pierwsze „koło” (strona 9 i następne)
  • Zadanie problemu z sumą arytmetyczną zadane w Math Meta
  • David Goldberg, „Co każdy informatyk powinien wiedzieć o arytmetyki zmiennoprzecinkowej”, ACM Computing Surveys 23 , 1 (1991-03), 5-48 doi> 10.1145 / 103162.103163 ( dostępna również wersja )
  • Przewodnik zmiennoprzecinkowy - co każdy programista powinien wiedzieć o arytmetyki zmiennoprzecinkowej
  • 0.30000000000000004.com porównuje arytmetykę zmiennoprzecinkową we wszystkich językach programowania
• Duplikat kanoniczny „zmiennoprzecinkowy jest niedokładny” (meta dyskusja na temat kanonicznej odpowiedzi na ten problem)
• Kilka pytań przepełnienia stosu, w tym
  • Dlaczego liczby dziesiętne nie mogą być przedstawione dokładnie w postaci binarnej?
  • Dlaczego liczby zmiennoprzecinkowe są niedokładne?
  • Czy matematyka zmiennoprzecinkowa jest zepsuta?
• Sztuczki matematyczne wyjaśnione przez Arthura T. Benjamina

. dplyr::near()jest kolejną opcją do testowania, czy dwa wektory liczb zmiennoprzecinkowych są równe.

Funkcja ma wbudowany parametr tolerancji: tol = .Machine$double.eps^0.5który można regulować. Domyślny parametr jest taki sam jak domyślny dla all.equal().