0

Could someone please explain why R has a problem with rounding/ceiling and how you can make sure that it calculates correctly?

N = 50e3
x = vector(length = N)
p = 0.995

ceiling(length(x)*(1-p))
251

ceiling(N*(1-p))
251

ceiling(N*0.005)
250

ceiling(50e3*0.005)
250
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • 3
    https://stackoverflow.com/questions/9508518/why-are-these-numbers-not-equal – rawr May 24 '21 at 17:22
  • I didn't think this was an exact duplicate because of the "and how can you make sure that it calculates correctly?" part ... – Ben Bolker May 24 '21 at 17:44

1 Answers1

3

@rawr links to Why are these numbers not equal? , which explains the fundamental issues with floating-point computation. In this case:

print((1-0.995) - 0.005)
## [1] 4.336809e-18

Because of this, (1-0.005)*5e4 is slightly greater than 250 (to see this you have to print((1-0.005)*5e4, digits=22), because R prints a rounded representation by default) so ceiling() pushes the answer up to 251.

In this particular case it looks like you can get the desired answer by rounding (1-p) to three decimal places, ceiling(N*round(1-p,3)) — but you should definitely read the linked answer and think about whether this solution will be robust for all of your needs.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453