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
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) +
