Вычислительные ошибки R, комбинации с большими числами создают неправильные числа только некоторые из времени - PullRequest
Новая
       6

Вычислительные ошибки R, комбинации с большими числами создают неправильные числа только некоторые из времени

0 голосов
/ 08 октября 2018

Я создал формулу комбинации в R для вычисления комбинаций больших чисел. 

  combination_1 <- function(n,r){
  n_0 <- n
  num <- 1
  denom <- factorial(r)
  for(i in 1:(r)){
    num <- num * n_0
    n_0 <- n_0-1
  }
  num/denom
}

Затем я применил эту формулу комбинации в циклическом тестировании, когда две величины равны, а затем распечатал соответствующиеЗначения n, m и r.

 both_large <- function(n,m,r){m*combination_1(n,r)}
 both_small <- function(n,m,r){factorial(m)*combination_1(n,(r-1))} 

 for(m in 3:8){
 for(n in 20:20000){
    for(r in 1:20){
      if(both_large(n,m,r) - both_small(n,m,r) == 0){
       cat('r = ', r, ', n = ', n,', m= ',m, '\n')
      }
    }
   }
}

Этот код, однако, работает только в некоторых случаях.Следующий вывод показан ниже с пропущенным значением при r = 9. 

r =  6 , n =  149 , m=  5 
r =  7 , n =  174 , m=  5 
r =  8 , n =  199 , m=  5 
r =  10 , n =  249 , m=  5 
r =  11 , n =  274 , m=  5

Определенно должно быть значение при r = 9, n = 224, m = 5;однако, когда я запускаю вычитание для этих конкретных значений, R вычисляет значение -2.Когда я запускаю его через Wolfram Alpha, он вычисляет значение 0. Я также нашел способ немного упростить мою формулу, и упрощенная версия также работает до 0.

Почему R не будетвычислить значение при r = 9, но затем правильно вычислить большие значения при r = 10 и r = 11?Это какая-то ошибка округления, и если это так, то почему она будет вычислять большие значения?Он также не будет вычислять другие значения.Это только первый случай, когда этого не происходит.

Спасибо!

1 Ответ

0 голосов
/ 09 октября 2018

Это проблема с двойной точностью .В статье мы видим, что двойные типы данных хранятся в 64 битах со следующей разбивкой:

  1. Бит знака: 1 бит
  2. Экспонент: 11 бит
  3. Значимость и точность: 53 бита (52 явно хранятся)

Преобразовав это в основание 10, мы увидим, что мы гарантированно получим не менее 15 десятичных цифр точности. 

log10(2^53 - 1)
[1] 15.95459

Мы можем увидеть, что это действие, наблюдая за странным поведением с простой арифметикой: 

options(scipen = 999)

1e16
[1] 10000000000000000

1e16 + 5
[1] 10000000000000004  ## incorrect.. should be 10000000000000005

На вашем примере r = 9, n = 224, m = 5 вместе с некоторыми print операторами внутри вашей функции combination_1 мы обнаружим виновника: 

combination_1_Verbose <- function(n,r,verbose = FALSE){
    n_0 <- n
    num <- 1
    denom <- factorial(r)
    for(i in 1:(r)){
        num <- num * n_0
        n_0 <- n_0-1
    }

    if (verbose) {
        print(num)
        print(log10(num))
    }

    num/denom
}

combination_1_Verbose(n, r - 1, TRUE)
[1] 5585745606995474432
[1] 18.74708
[1] 138535357316356

Мы выполняем арифметику из более чем 18 цифр ... вне пределов нашей точности.

Что также не очевидно, так это тот факт, что возвращаемое значение не совсем 138535357316356,Используя аргумент digits для print, мы фактически видим, что возвращаемое значение не является целым числом. 

print(combination_1_Verbose(n, r - 1), digits = 22)
[1] 138535357316356.015625

Это в конечном итоге является источником вашей ошибки.Если мы возьмем .015625 и умножим на factorial(m) = 120, получим: 

.015625 * 120
[1] 1.875

Это округляется до 2, что является разницей в нашем чеке.

Мы можемисправляет это поведение, используя библиотеку с множественной точностью gmp: 

library(gmp)
combination_1_GMP <- function(n,r,verbose = FALSE){
    n_0 <- as.bigz(n)
    num <- as.bigz(1)
    denom <- factorialZ(r)
    for(i in 1:(r)){
        num <- mul.bigz(num, n_0)
        n_0 <- sub.bigz(n_0, 1)
    }

    if (verbose) {
        print(num)
        print(log10(num))
    }

    as.bigz(num/denom)
}

combination_1_GMP(n, r-1, TRUE)
Big Integer ('bigz') :
[1] 5585745606995473920
[1] 18.74708
Big Integer ('bigz') :
[1] 138535357316356

В исходной функции num было 5585745606995474432, тогда как в нашем примере gmp мы получили 5585745606995473920.Обратите внимание, что разница меньше 500, число из 3 цифр.Это имеет смысл, поскольку наш номер превышает 18 цифр, и, как мы уже говорили выше, нам гарантируется только 15 общих цифр точности (т. Е. 18 - 3 = 15).

В качестве альтернативы, мы могли бы round конечный результат.Я не рекомендовал бы эту опцию, если точность абсолютно необходима, поскольку все еще есть некоторые значения n, m, r, которые все еще будут зависеть от двойной точности.Это работает в этом примере, хотя: 

combination_1_Round <- function(n,r){
    n_0 <- n
    num <- 1
    denom <- factorial(r)
    for(i in 1:(r)){
        num <- num * n_0
        n_0 <- n_0-1
    }
    round(num/denom)
}

both_large_r <- function(n,m,r){m*combination_1_Round(n,r)}
both_small_r <- function(n,m,r){factorial(m)*combination_1_Round(n,(r-1))}

both_large_r(n,m,r) - both_small_r(n,m,r)
[1] 0

И, наконец, ваш лучший вариант - переписать ваш алгоритм, чтобы сохранить числа в пределах двойной точности. 

combination_1_Improved <- function(n,r){

    denom <- num <- 1
    i <- (n - r + 1)

    for (denom in 1:r) {
        num <- num * i;
        num <- num / denom;
        i <- i + 1
    }

    num
}

print(combination_1_Improved(n,r-1), digits = 22)
[1] 138535357316356
...