From: Charles C. Berry <cberry_at_tajo.ucsd.edu>

Date: Fri, 14 Sep 2007 10:18:06 -0700

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri 14 Sep 2007 - 17:22:48 GMT

Date: Fri, 14 Sep 2007 10:18:06 -0700

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 Diegohttp://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri 14 Sep 2007 - 17:22:48 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Fri 14 Sep 2007 - 18:40:52 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*