From a question on Stackoverflow: http://stackoverflow.com/questions/23365599/how-to-determine-the-standard-normal-loss-function-lz/

2014-04-30 , Martin Maechler, ETH Zurich

Note that I define L(x) already amended to provide the limit \( L(\infty) = 0 \), as well:

L <- function(x) ifelse(x == Inf, 0, dnorm(x) - x * pnorm(x, lower.tail = FALSE))

curve(L, -3, 3)
abline(h = 0, v = 0, lty = 3, col = adjustcolor("gray", 0.8))  # looks fine
points(0, 1/sqrt(2 * pi), col = adjustcolor(2, 0.5), pch = 7)

curve(L, -1, 10)

curve(L, -1, 10, log = "y")

curve(L, -10, 50, log = "y", yaxt = "n")

## Warning: 20 y values <= 0 omitted from logarithmic plot

sfsmisc::eaxis(2)

Now do the y-axis log transformed: partly want to use to find good interval for root finding in inversion:

curve(log(L(x)), -10, 50)

The original simple proposition

Linv.v1 <- function(y) uniroot(function(x) L(x) - y, interval = c(-10, 10))$root

works in typical cases, but fails for vector input or for very small y or y values larger than 10.

Here, inspired by the above plots, I do the inversion fast on the log-scale when “on the right”, and care for the extreme cases:

.Linv <- function(y, tol = 1e-05, ...) {
    stopifnot(length(y) == 1, is.numeric(y))
    if (y <= 0) 
        return(Inf)
    if (y < 1/sqrt(2 * pi)) {
        ## result > 0, work on log scale
        l.y <- log(y)
        if (l.y < -744.440071921381) 
            return(39)  # 38.568041815493341312
        uniroot(function(x) log(L(x)) - l.y, interval = c(0, 1.55 * sqrt(-l.y)), 
            tol = tol, ...)$root
    } else {
        ## y >= 0.3989 ==> result <= 0
        uniroot(function(x) L(x) - y, interval = c(min(-1, -1.01 * y), 0), tol = tol, 
            ...)$root
    }
}

.Linv(2)  # ok

## [1] -1.991

Now the vectorizing version : (you need, e.g., for curve(.):

Linv <- function(x, ...) {
    r <- -x  # numerically exact for x >= 8;  use .Linv(.) for x < 8 :
    ii <- x < 8
    r[ii] <- vapply(x[ii], .Linv, 1)
    r
}
curve(Linv, -1, 10)

looks ok, but I know I can do better, notably plotting *both \( L \) and its inverse \( L^{-1}() = \) Linv in the same plot and showing that they are mirror images of each other:

p.L.and.inv <- function(from, to, n = 2049) {
    c1 <- "black"
    c2 <- "tomato"
    xy <- curve(Linv, from, to, n = n, asp = 1, col = c1, lwd = 2, ylab = "", 
        main = quote(list(L(x) ~ "and" ~ {
            L^-1
        }(x), ~~~L(x) == phi(x) - x * (1 - Phi(x)))))
    curve(L, add = TRUE, col = c2, lwd = 2, n = n)
    abline(h = 0, v = 0, lty = 2, col = adjustcolor("gray", 0.8))
    abline(0, 1, col = "dark green", lty = 3)
    abline(0, -1, col = "blue3", lty = 4)
    legend("top", expression(L(x), {
        L^-1
    }(x)), lwd = 2, col = c(c2, c1), bty = "n")
    invisible(xy)
}

and now a first plot, and a zoom in of it:

p.L.and.inv(-12, 12)

p.L.and.inv(-0.5, 2)

Well, now that ended up in a little lesson about functions and inverses, using R.

_______That's all, folks!_______