I'd compute this in the log-scale (taking also advantage of the 'log' and 'log.p' arguments of dnorm() and pnorm(), respectively), and then transform back, e.g.,
fn1 <- function(B){ -(pnorm(B) * dnorm(B) * B + dnorm(B)^2)/pnorm(B)^2 } fn2 <- function(B){ p1 <- dnorm(B, log = TRUE) + log(-B) - pnorm(B, log.p = TRUE) p2 <- 2 * (dnorm(B, log = TRUE) - pnorm(B, log.p = TRUE)) exp(p1) - exp(p2) } fn1(c(-15, -25, -35, -55, -105)) fn2(c(-15, -25, -35, -55, -105)) I hope it helps. Best, Dimitris ---- Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm Quoting Ott Toomet <[EMAIL PROTECTED]>: > 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, > -1) > > 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, > Ott > > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > > Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html