Re: [Rd] Reciprocal Mill's Ratio

From: Charles C. Berry <cberry_at_tajo.ucsd.edu>
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, 2232.
>
> 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



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.