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.hard...@wlandres.net> Date: 19-Jun-2012 Time: 00:07:11 This message was sent by XFMail ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
