[Rd] Accuracy (PR#14139)

[savicky at cs.cas.cz]

On Mon, Dec 14, 2009 at 06:10:16PM +0100, bersch at lycos.com wrote:
> 
> > pnorm(1.35,0,1)
> [1] 0.911492
> > pnorm(1.36,0,1)
> [1] 0.913085
> 
> > options(digits=4)
> 
> > pnorm(1.35,0,1)
> [1] 0.9115
> > pnorm(1.36,0,1)
> [1] 0.913          rounding error?

The technical explanation is as follows. If options(digits=k) is set,
then the number of significant digits for printing a single number x is 
determined as min(k, d), where d is the minimum number of digits,
for which the relative error of the printed number is less than 10^-k.

If we have
  x <- 0.913085
  y <- 0.913
then the relative error of y as an approximation of x is
  abs(y - x)/x # [1] 9.3091e-05
Since this is less than 10^-4, the 3 digit precision is chosen
for printing x.

A safer way of rounding is to use functions round() and signif().
For example,
  round(x, digits=4) # [1] 0.9131

I do not know the history of the R printing algorithm. It is designed
primarily for printing vectors, where the rules are more complicated
to achieve a good unified format for all numbers. May be, someone else
can say more about it. The above analysis may be obtained by inspecting
the R source code.

Petr Savicky.