Creating Mills Ratio in R for large values

Question

I'm using R to create a function, that amongst others uses Mills Ratio (See here). This is not a complicated formula, and at first I just programmed it like this:

mill <- function(x) {
  return((1 - pnorm(x)) / dnorm(x))
}

I soon found out however, that for very large values (x >= 9) of x , this function returns zero. Even more dramatic, at around x >= 37, it starts returning NaN , which really messes up my stuff.

Following the article, for now I've changed the function into this:

mill <- function(x) {
  if (x >= 9) {
    return(1 / x)
  } else {
    return((1 - pnorm(x)) / dnorm(x))
  }
}

This seems to work. I use this function to calculate a vector however, and when I use simulation to find the same vector, I get more or less the same answer, only it's always a bit off..

I think this has to do with my implementation of Mills Ratio, since the rest of the function is just exponentials, which R should have no trouble with.

I want to ask you guys if there is any way to solve this problem: to either implement this function better, or give me another way to find the Mills Ratio (perhaps through integration of some sorts, but wouldn't I run into the same issues there?). Thank you kindly for any help you can provide!

Answer 1

I would make two changes to your original mill function.

1. Change 1-pnorm(x) to pnorm(lower.tail=FALSE)
2. Use log's and take exponentials if needed.

So this gives

new_mill = function(x) 
    pnorm(x, lower.tail=FALSE, log.p=TRUE) - dnorm(x, log=TRUE)

So

R> exp(new_mill(10))
[1] 0.09903
R> exp(new_mill(40))
[1] 0.02498

Using a plot as a sanity check

x = seq(0, 10, 0.001)
plot(x, exp(new_mill(x)), type="l")
lines(x, mill(x), col=2)

gives