[R] workaround for numeric problems

Ott Toomet otoomet at ut.ee
Sun Jul 2 12:00:10 CEST 2006

Dear R-people,

I have to compute 

C - -(pnorm(B)*dnorm(B)*B + dnorm(B)^2)/pnorm(B)^2

This expression seems to be converging to -1 if B approaches to -Inf
(although I am unable to prove it).  R has no problems until B equals
around -28 or less, where both numerator and denominator go to 0 and
you get NaN. A simple workaround I did was

C <- ifelse(B > -25,
           -(pnorm(B)*dnorm(B)*B + dnorm(B)^2)/pnorm(B)^2,

It works well for me (32bit intel/linux platform).  But what about
other processors/platforms/compilator options?  Are there any better
ways for finding out at which values the numerical problems start?
Can one derive something from .Machine$double.eps (but what about the
precison of dnorm and other analytic functions)?

Thanks in advance,

More information about the R-help mailing list