Apa cara yang benar / standar untuk memeriksa apakah perbedaan lebih kecil dari presisi mesin?

Saya sering berakhir dalam situasi di mana perlu untuk memeriksa apakah perbedaan yang diperoleh di atas presisi mesin. Sepertinya untuk tujuan ini R memiliki variabel berguna: .Machine$double.eps. Namun ketika saya beralih ke kode sumber R untuk panduan tentang menggunakan nilai ini saya melihat beberapa pola yang berbeda.

Contohnya

Berikut beberapa contoh dari statsperpustakaan:

t.test.R

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

chisq.test.R

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

mengintegrasikan

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

dll.

Pertanyaan

1. Bagaimana seseorang dapat memahami alasan di balik semua berbeda 10 *, 100 *, 50 *dan sqrt()pengubah?
2. Apakah ada pedoman tentang cara .Machine$double.epsmenyesuaikan perbedaan karena masalah presisi?

Presisi mesin doubletergantung pada nilai saat ini. .Machine$double.epsmemberikan presisi saat nilainya 1. Anda dapat menggunakan fungsi C nextAfteruntuk mendapatkan presisi mesin untuk nilai lainnya.

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

Menambahkan nilai ake nilai btidak akan berubah bketika aadalah <= setengah dari itu mesin presisi. Memeriksa apakah perbedaannya lebih kecil daripada presisi mesin dilakukan <. Pengubah mungkin mempertimbangkan kasus-kasus tertentu seberapa sering penambahan tidak menunjukkan perubahan.

Dalam R , presisi mesin dapat diperkirakan dengan:

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

Setiap doublenilai mewakili rentang. Untuk tambahan sederhana, kisaran hasil tergantung pada pengubahan setiap musim dan juga jumlah dari jumlah mereka.

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

Untuk precisi yang lebih tinggi Rmpfrbisa digunakan.

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

Dalam hal ini dapat dikonversi ke integer gmpdapat digunakan (apa yang ada di Rmpfr).

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

Definisi machine.eps: ini adalah nilai terendah  eps yang  1+eps bukan 1

Sebagai aturan praktis (dengan asumsi representasi floating point dengan basis 2):
Ini epsmembuat perbedaan untuk rentang 1 .. 2,
untuk rentang 2 .. 4 presisinya 2*eps
dan seterusnya.

Sayangnya, tidak ada aturan praktis yang baik di sini. Ini sepenuhnya ditentukan oleh kebutuhan program Anda.

Di R kita memiliki semua. Sama dengan cara bawaan untuk menguji perkiraan persamaan. Jadi Anda bisa menggunakan mungkin sesuatu seperti (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 memiliki sejumlah pencocokan floating point untuk perbandingan presisi ganda, termasuk DoubleEqdan DoubleNear. Anda dapat menggunakannya dalam pencocokan array seperti ini:

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

Memperbarui:

Numerical Recipes memberikan derivasi untuk menunjukkan bahwa dengan menggunakan selisih selisih satu sisi, sqrtadalah pilihan langkah-ukuran yang baik untuk perkiraan perbedaan turunan yang terbatas.

Situs artikel Wikipedia, Numerical Recipes, edisi ke-3, Bagian 5.7, yang merupakan halaman 229-230 (sejumlah terbatas tampilan halaman tersedia di http://www.nrbook.com/empanel/ ).

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

Aritmatika titik apung IEEE ini adalah batasan aritmatika komputer yang terkenal dan dibahas di beberapa tempat:

• FAQ di R memiliki seluruh pertanyaan yang didedikasikan untuknya: R FAQ 7.31
  • The R Inferno oleh Patrick Burns mempersembahkan "Lingkaran" pertama untuk masalah ini (halaman 9 dan seterusnya)
  • Aritmatika Jumlah Bukti yang diminta dalam Matematika Meta
  • David Goldberg, "Apa Yang Harus Diketahui Setiap Ilmuwan Mengenai Aritmatika Poin-Mengapung," ACM Computing Survey 23 , 1 (1991-03), 5-48 doi> 10.1145 / 103162.103163 ( revisi juga tersedia )
  • The Floating-Point Guide - Yang Harus Diketahui Setiap Programmer Tentang Aritmatika Floating-Point
  • 0.30000000000000004.com membandingkan aritmatika titik apung di seluruh bahasa pemrograman
• Duplikat kanonik untuk "floating point isaccurate" (diskusi meta tentang jawaban kanonik untuk masalah ini)
• Beberapa pertanyaan Stack Overflow termasuk
  • Mengapa angka desimal tidak dapat direpresentasikan secara tepat dalam biner?
  • Mengapa angka floating point tidak akurat?
  • Apakah matematika floating point rusak?
• Trik Matematika Dijelaskan oleh Arthur T. Benjamin

. dplyr::near()adalah pilihan lain untuk menguji apakah dua vektor angka floating point sama.

Fungsi ini memiliki parameter toleransi bawaan: tol = .Machine$double.eps^0.5yang dapat disesuaikan. Parameter default sama dengan default untuk all.equal().