Wie kann überprüft werden, ob der Unterschied kleiner als die Maschinengenauigkeit ist?

36

Ich lande oft in Situationen, in denen überprüft werden muss, ob der ermittelte Unterschied über der Maschinengenauigkeit liegt. Scheint für diesen Zweck R hat eine handliche Variable:.Machine$double.eps . Wenn ich mich jedoch an den R-Quellcode wende, um Richtlinien zur Verwendung dieses Werts zu erhalten, werden mehrere unterschiedliche Muster angezeigt.

Beispiele

Hier einige Beispiele aus stats Bibliothek:

t.test.R

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

chisq.test.R

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

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

usw.

Fragen

  1. Wie kann man die Gründe für all diese verschiedenen verstehen 10 *, 100 *, 50 *und sqrt()Modifikatoren?
  2. Gibt es Richtlinien .Machine$double.epszur Anpassung von Unterschieden aufgrund von Präzisionsproblemen?
Karolis Koncevičius
quelle
6
Beide Beiträge kommen daher zu dem Schluss, dass "das angemessene Maß an Sicherheit" von Ihrer Bewerbung abhängt. Als Fallstudie können Sie diesen Beitrag auf R-devel überprüfen . "Aha! 100-fache Maschinenpräzision in nicht allzu viel, wenn die Zahlen selbst zweistellig sind." (Peter Dalgaard, Mitglied des R Core-Teams)
Henrik
1
@ KarolisKoncevičius, ich denke nicht, dass es so einfach ist. Dies hängt mit den allgemeinen Fehlern in der Gleitkomma-Mathematik und der Anzahl der Operationen zusammen, die Sie mit ihnen ausführen. Wenn Sie einfach mit Gleitkommazahlen vergleichen, verwenden Sie double.eps. Wenn Sie mehrere Operationen mit einer Gleitkommazahl ausführen, sollte sich auch Ihre Fehlertoleranz anpassen. Deshalb gibt Ihnen all.equal ein toleranceArgument.
Joseph Wood
1
Schauen Sie sich auch die Implementierung der nextafter-Funktionalität in R an , um die nächstgrößere doppelte Zahl zu erhalten.
GKi

Antworten:

4

Die Maschinengenauigkeit für doublehängt von ihrem aktuellen Wert ab. .Machine$double.epsGibt die Genauigkeit an, wenn die Werte 1 sind. Mit der C-Funktion können nextAfterSie die Maschinengenauigkeit für andere Werte ermitteln.

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

Das Hinzufügen von Wert azu Wert bändert sich nicht, bwenn adie <= Hälfte der Maschinengenauigkeit erreicht ist. Es wird geprüft, ob der Unterschied geringer ist als die Maschinengenauigkeit< . Die Modifikatoren können typische Fälle berücksichtigen, in denen eine Addition keine Änderung zeigte.

In R kann die Maschinengenauigkeit geschätzt werden mit:

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

Jeder doubleWert repräsentiert einen Bereich. Für eine einfache Addition hängt der Bereich des Ergebnisses von der Neuausrichtung jedes Summanden und auch vom Bereich ihrer Summe ab.

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

Für eine höhere Präzision Rmpfrkönnte verwendet werden.

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

Falls es in eine Ganzzahl konvertiert werden gmpkönnte, könnte es verwendet werden (was in Rmpfr ist).

library(gmp)
as.bigz("2")^54 * 10 - 29 + 41 - (as.bigz("2")^54 * 10 + 59)
#[1] -47
GKi
quelle
Vielen Dank. Ich denke, das ist eine viel bessere Antwort. Es zeigt viele Punkte gut. Das einzige, was mir noch etwas unklar ist, ist: Kann man sich die Modifikatoren (wie * 9 usw.) selbst einfallen lassen? Und wenn ja wie ...
Karolis Koncevičius
Ich denke, dieser Modifikator entspricht dem Signifikanzniveau in der Statistik und erhöht sich um die Anzahl der Operationen, die Sie in Kombination mit dem gewählten Risiko durchgeführt haben, um einen korrekten Vergleich abzulehnen.
GKi
3

Definition einer machine.eps: Dies ist der niedrigste Wert,   für den dies   nicht der Fall ist eps1+eps1

Als Faustregel (unter der Annahme einer Gleitkommadarstellung mit Basis 2):
Dies epsmacht den Unterschied für den Bereich 1 .. 2,
für den Bereich 2 .. 4 ist die Genauigkeit2*eps
und so weiter.

Leider gibt es hier keine gute Faustregel. Es wird ganz von den Bedürfnissen Ihres Programms bestimmt.

In R haben wir all.equal als eingebauten Weg, um die ungefähre Gleichheit zu testen. Also könntest du vielleicht so etwas gebrauchen (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 verfügt über eine Reihe von Gleitkomma-Matchern für Vergleiche mit doppelter Genauigkeit, einschließlich DoubleEqund DoubleNear. Sie können sie in einem Array-Matcher wie folgt verwenden:

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

Aktualisieren:

Numerische Rezepte liefern eine Ableitung, um zu demonstrieren, dass unter Verwendung eines einseitigen Differenzquotienten sqrt eine gute Wahl der Schrittgröße für endliche Differenznäherungen von Derivaten ist.

Die Wikipedia-Artikelseite Numerical Recipes, 3. Auflage, Abschnitt 5.7, Seiten 229-230 (eine begrenzte Anzahl von Seitenaufrufen ist unter http://www.nrbook.com/empanel/ verfügbar ).

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

Diese IEEE-Gleitkomma-Arithmetik ist eine bekannte Einschränkung der Computer-Arithmetik und wird an mehreren Stellen diskutiert:

. dplyr::near()ist eine weitere Option zum Testen, ob zwei Vektoren von Gleitkommazahlen gleich sind.

Die Funktion verfügt über einen eingebauten Toleranzparameter tol = .Machine$double.eps^0.5, der angepasst werden kann. Der Standardparameter ist der gleiche wie der Standardparameter für all.equal().

Sreeram Nair
quelle
2
Danke für die Antwort. Im Moment denke ich, dass dies zu minimal ist, um eine akzeptierte Antwort zu sein. Es scheint nicht die beiden Hauptfragen aus dem Beitrag zu beantworten. Zum Beispiel heißt es "es wird durch die Bedürfnisse Ihres Programms bestimmt". Es wäre schön, ein oder zwei Beispiele für diese Aussage zu zeigen - vielleicht ein kleines Programm und wie Toleranz dadurch bestimmt werden kann. Möglicherweise mit einem der genannten R-Skripte. Hat all.equal()auch seine eigene Annahme als Standardtoleranz gibt es sqrt(double.eps)- warum ist es der Standard? Ist es eine gute Faustregel sqrt()?
Karolis Koncevičius
Hier ist der Code, mit dem R das EPS berechnet (in ein eigenes Programm extrahiert). Außerdem habe ich die Antwort mit zahlreichen Diskussionspunkten aktualisiert, die ich zuvor durchlaufen hatte. Hoffe das gleiche hilft dir besser zu verstehen.
Sreeram Nair
Eine aufrichtige +1 für alle Anstrengungen. Aber zum gegenwärtigen Zeitpunkt kann ich die Antwort immer noch nicht akzeptieren. Es scheint ein bisschen unerreichbar mit vielen Referenzen zu sein, aber in Bezug auf die tatsächliche Antwort auf die 2 gestellten Fragen: 1) wie man 100x, 50x usw. Modifikatoren in R- stats::Quelle versteht und 2) was sind die Richtlinien; Die Antwort ist ziemlich dünn. Der einzig zutreffende Satz scheint der Verweis aus "Numerical Recipes" zu sein, wonach sqrt () ein guter Standard ist, was meiner Meinung nach wirklich zutreffend ist. Oder vielleicht fehlt mir hier etwas.
Karolis Koncevičius