¿Cuál es la forma correcta / estándar de verificar si la diferencia es menor que la precisión de la máquina?

A menudo termino en situaciones en las que es necesario verificar si la diferencia obtenida está por encima de la precisión de la máquina. Parece que para este propósito R tiene una variable útil: .Machine$double.eps. Sin embargo, cuando recurro al código fuente R para obtener pautas sobre el uso de este valor, veo múltiples patrones diferentes.

Ejemplos

Aquí hay algunos ejemplos de la 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.

Preguntas

1. ¿Cómo entender el razonamiento detrás de todas esas diferentes 10 *, 100 *, 50 *y sqrt()modificadores?
2. ¿Existen pautas sobre el uso .Machine$double.epspara ajustar las diferencias debido a problemas de precisión?

La precisión de la máquina doubledepende de su valor actual. .Machine$double.epsproporciona la precisión cuando los valores son 1. Puede usar la función C nextAfterpara obtener la precisión de la máquina para otros 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

Agregar valor aa valor bno cambiará bcuando asea ​​la <= mitad de la precisión de su máquina. Verificando si la diferencia es más pequeña que con la precisión de la máquina <. Los modificadores podrían considerar casos típicos con qué frecuencia una adición no mostró un cambio.

En R, la precisión de la máquina se puede estimar con:

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 representa un rango. Para una simple suma, el rango del resultado depende del reange de cada sumando y también del rango de su suma.

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 mayor precisión se Rmpfrpodría utilizar.

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

En caso de que pudiera convertirse a entero, gmppodría usarse (lo que está en Rmpfr).

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

Definición de una máquina.eps: es el valor más bajo  eps para el que  1+eps no 1

Como regla general (suponiendo una representación de punto flotante con base 2):
Esto epshace la diferencia para el rango 1 .. 2,
para el rango 2 .. 4 la precisión es 2*eps
y así sucesivamente.

Desafortunadamente, no hay una buena regla general aquí. Está completamente determinado por las necesidades de su programa.

En R tenemos all.equal como una forma integrada para probar la igualdad aproximada. Entonces podrías 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

Google simulacro tiene una serie de comparadores de punto flotante para comparaciones de doble precisión, incluyendo DoubleEqy DoubleNear. Puede usarlos en un comparador de matrices como este:

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

Actualizar:

Las recetas numéricas proporcionan una derivación para demostrar que usar un cociente de diferencia unilateral sqrtes una buena opción de tamaño de paso para aproximaciones de derivadas de diferencias finitas.

El sitio del artículo de Wikipedia Numerical Recipes, 3ra edición, Sección 5.7, que está en las páginas 229-230 (hay un número limitado de visitas a la página en http://www.nrbook.com/empanel/ ).

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

Esta aritmética de coma flotante IEEE es una limitación bien conocida de la aritmética informática y se discute en varios lugares:

• Las preguntas frecuentes en R tienen toda una pregunta dedicada: preguntas frecuentes R 7.31
  • El Infierno R de Patrick Burns dedicó el primer "Círculo" a este problema (página 9 en adelante)
  • Problema de prueba de suma aritmética en Meta Meta
  • David Goldberg, "Lo que todo informático debe saber sobre la aritmética de punto flotante", ACM Computing Surveys 23 , 1 (1991-03), 5-48 doi> 10.1145 / 103162.103163 ( revisión también disponible )
  • La guía de punto flotante: lo que todo programador debe saber sobre la aritmética de punto flotante
  • 0.30000000000000004.com compara la aritmética de coma flotante en los lenguajes de programación
• Duplicado canónico para "punto flotante es inexacto" (una meta discusión sobre una respuesta canónica para este problema)
• Varias preguntas de desbordamiento de pila que incluyen
  • ¿Por qué los números decimales no se pueden representar exactamente en binario?
  • ¿Por qué los números de coma flotante son inexactos?
  • ¿Se han roto las matemáticas de coma flotante?
• Trucos matemáticos explicados por Arthur T. Benjamin

. dplyr::near()es otra opción para probar si dos vectores de números de coma flotante son iguales.

La función tiene un parámetro de tolerancia incorporado: tol = .Machine$double.eps^0.5que se puede ajustar. El parámetro predeterminado es el mismo que el predeterminado para all.equal().