ley - ¿Por qué R realiza incorrectamente la suma aquí?

didi uber

Question

(Si la respuesta es, de hecho, cómo R usa números binarios para almacenar pequeños decimales, aún así encantaría escuchar detalles sobre por qué esto sucede a alguien con conocimiento sobre esto)

Estoy usando R para calcular una suma, específicamente esta:

Aquí está mi código R:

r = 21 #Number of times the loop executes 
n = 52 
sum = 0 #Running total 
for(k in 0:r){
  sum = sum + (1-(choose(2^(k), n)*(factorial(n)/(2^(n*k)))))
  print(sum)
}

Si miras la salida, notarás:

[1] 1
[1] 2
...
[1] 11.71419
[1] 11.71923
[1] 11.72176
[1] 12.72176
[1] 13.72176

¿Por qué comienza a incrementarse en 1 después de la decimonovena iteración?

¿Hay otros motores computacionales disponibles libremente que serían más adecuados para esta tarea?

Answer 1

Todos los cálculos están usando números de coma flotante aquí, que generalmente no es muy adecuado para factores y potencias de elevación (porque tales valores se hacen rápidamente muy grandes, lo que hace que el cálculo sea menos preciso).

Por favor, inténtalo:

> factorial(52)/(2^(52*19))
[1] 3.083278e-230
> factorial(52)/(2^(52*20))
[1] 0

y compara con Pari-GP:

? /p 256
realprecision = 269 significant digits (256 digits displayed)
? precision( (52!) / 2^(52*19) + .  , 24)
%1 = 3.08327794368108826958435659488643289724 E-230
? precision( (52!) / 2^(52*20) + .  , 24)
%2 = 6.8462523287873017654100595727727116496 E-246

Answer 2

Si está buscando una forma de evitar el problema de desbordamiento / desbordamiento que tiene, puede usar registros (para mantener las magnitudes razonables para los cálculos intermedios) y luego exponer al final:

options(digits=20)

for(k in 0:r){
  sum = sum + (1 - exp(lchoose(2^k, n) + log(factorial(n)) - k*n*log(2)))
  print(paste0(k,": ",sum))
}

[1] "0: 1"
[1] "1: 2"
[1] "2: 3"
...
[1] "19: 11.7217600143979"
[1] "20: 11.7230238079842"
[1] "21: 11.7236558993777"

Para verificar que esto sea correcto, ejecuté la suma original (sin tomar registros) en Mathematica y obtuve los mismos resultados con 12 decimales.

Aunque puede resolver el problema en R , si quiere usar un sistema de álgebra computacional (que le permite hacer cálculos simbólicos y cómputos exactos), Sage es de código libre y de código abierto.