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. 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.
Dear Murray, I think you might have to update R in first instance and provide a reproducible example in second place so that people might help you. Regards, Simone On 13/09/2007, maj@stats.waikato.ac.nz <maj@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. > > 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@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel >-- ______________________________________________________ Simone Giannerini Dipartimento di Scienze Statistiche "Paolo Fortunati" Universita' di Bologna Via delle belle arti 41 - 40126 Bologna, ITALY Tel: +39 051 2098262 Fax: +39 051 232153 http://www2.stat.unibo.it/giannerini/ ______________________________________________________ [[alternative HTML version deleted]]
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