[R] (1-1e-100)==1 true?

(Ted Harding) Ted.Harding at wlandres.net
Tue Jun 19 01:07:15 CEST 2012


On 18-Jun-2012 21:26:41 whf1984911 wrote:
> Hi, 
> 
> This problems has bothered me for the lase couple of hours.  
> 
>> 1e-100==0
> [1] FALSE
>> (1-1e-100)==1
> [1] TRUE
> 
> How can I tell R that 1-1e-100 does not equal to 1,  actually,
> I found out that
>  > (1-1e-16)==1
> [1] FALSE
>> (1-1e-17)==1
> [1] TRUE
> 
> The reason I care about this is that I was try to use qnorm()
> in my code, for example,
> 
>> qnorm(1e-100)
> [1] -21.27345
> 
> and if I want to find qnorm(x) where x is very close to 1, say
> x=1-1e-100, then you would think using qnorm(1-x, lower.tail=F)
> would give me something other than INF, but that does not work
> since R would recognize x==1 in this case and therefore, 1-x==0,
> so qnorm(1-x, lower.tail=F) will give me INF which is what I try
> to avoid in my code.
> 
> How could get around this, any suggestions?
> 
> Thanks,
> Haifeng Wu
> Graduate Student
> University of South Carolina

You are working close to, and also beyond, the boundary of what
R can internally represent, so anomalies like this will frequently
turn up.

For the sort of extreme values which you are using, it is better
to make use of asymptotic formulae, and not to use R's standard
functions such as qnorm, pnorm, ...

For example, there is Mills Ratio, for which (in RF notation)

  x*(1 - pnorm(x))/dnorm(x) --> 1 as x --> Inf

So, for large x, proceed as though it was exactly 1.

Then

  log(x) + log(1 - pnorm(x)) = log(dnorm(x))
  = -0.5*log(2*pi) - 0.5*x^2

so, with (e.g.) pnorm(x) = 1-1e-100, you have

  log(x) + 0.5*x^2 = 100 - 0.5*log(2*pi)

More generally, with pnorm(x) = 1-1e-X where X is large,

  log(x) + 0.5*x^2 = X - 0.5*log(2*pi)

so you now need to write a function, say Qnorm(X), which
solves this equation for x, given X.

The above is only a suggestion. Probably others know a better
approach!

Ted.

-------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at wlandres.net>
Date: 19-Jun-2012  Time: 00:07:11
This message was sent by XFMail



More information about the R-help mailing list