[Rd] Reciprocal Mill's Ratio

[Charles C. Berry]

On Fri, 14 Sep 2007, maj at stats.waikato.ac.nz wrote:

> I believe that this may be more appropriate here in r-devel than in r-help.
>
> The normal hazard function, or reciprocal Mill's Ratio, may be obtained
> in R as dnorm(z)/(1 - pnorm(z)) or, better, as dnorm(z)/pnorm(-z) for
> small values of z. The latter formula breaks dowm numerically for me
> (running R 2.4.1 under Windows XP 5.1 SP 2) for values of z near 37.4
> or greater.


OK,

> mr <- function(z )dnorm( z )/ ( pnorm(z,lower=FALSE) )
> mr.2 <- function(z) exp(dnorm( z, log=TRUE ) - (pnorm(z, lower=FALSE, log=TRUE )))
>
> mr(seq(10,100,by=10)) # breaks before 40
  [1] 10.09809 20.04975 30.03326      NaN      NaN      NaN      NaN      NaN      NaN      NaN
>
> mr.2(seq(10,100,by=10)) #seems robust
  [1]  10.09809  20.04975  30.03326  40.02497  50.01998  60.01666  70.01428  80.01250  90.01111 100.01000
> 
>
> mr.2(1e4)
[1] 10000
> mr.2(1e5)
[1] 99999.97
> mr.2(1e6)
[1] 999980.2
> mr.2(1e7)
[1] 9990923
> mr.2(1e8) # Oh well!
[1] 65659969
>

I guess you get large than 1e4, you should just use the asymptotic result.

HTH,

Chuck


>
> Looking at the pnorm documentation I see that it is based on Cody (1993)
> and thence, going one step further back, on Cody (1969). Consulting
> Cody (1969) I see that the algorithm for pnorm(z) [or actually erf(z)]
> is actually based on rational function approximations for the
> reciprocal Mill's ratio itself, as I rather expected.
>
> I wonder if anyone has dug out a function for the reciprocal Mill's
> ratio out of the pnorm() code? Anticipating the obvious response I don't
> believe that this would be one of the things I might be good at!
>
> Murray Jorgensen
>
> References
>
> Cody, W. D. (1993)
> Algorithm 715: SPECFUN – A portable FORTRAN package of special function
> routines and test drivers.
> ACM Transactions on Mathematical Software 19, 22–32.
>
> Cody, W. D. (1969)
> Rational Chebyshev Approximations for the Error Function
> Mathematics of Computation, Vol. 23, No. 107. (Jul., 1969), pp. 631-637.
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901