[Rd] Reciprocal Mill's Ratio
Charles C. Berry
cberry at tajo.ucsd.edu
Fri Sep 14 19:18:06 CEST 2007
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
More information about the R-devel
mailing list