Laorie - 1 year ago 99
R Question

# Standard Normal Quantile Function Integration in R

I need to compute a division of integrals, where the function q_alpha(z) is the quantile function of a standard normal distribution.

I got a question regarding the denominator. As the normal standard distribution has Homoscedasticity, it is simmetric, continuous, etc.The integration of the denominator term its simple? I just need to elevated to the square each quantile of this function and proceed to the calculation? Right?

This is my code in R:

``````library(Bolstad)

thau=1:99/100
z.standard.quantile=qnorm(thau,0,1)
z.standard.quantile.square=qnorm(thau,0,1)^2

sintegral(thau[1:50],z.standard.quantile[1:50])\$value/sintegral(thau[1:50], z.standard.quantile.square[1:50])\$value
``````

The result is:
`-0.8676396`

There is no problem in taking the square of `qnorm`, but `qnorm` is unbounded on `[0, 0.5]` (note `qnorm(0)` is `-Inf`) so the integral is not finite.

My second thought is that there is actually no need to use `Bolstad::sintegral` (Simpson's rule); the R base function `integrate` is sufficient. Or, we can discretize `qnorm` and use Trapezoidal rule because `qnorm` is a smooth function which can be well approximated by linear interpolation.

I will write a function evaluating the ratio of integral in your question, but lower bounded on `l`:

``````## using `integrate`
f1 <- function (l) {
a <- integrate(qnorm, lower = l, upper = 0.5)\$value
b <- integrate(function (x) qnorm(x) ^ 2, lower = l, upper = 0.5)\$value
a / b
}

## using Trapezoidal rule, with `n` division on interval `[l, 0.5]`
f2 <- function (l, n) {
x <- seq(l, 0.5, length = n)
delta <- x[2] - x[1]
y1 <- qnorm(x)
y2 <- y1 ^ 2
a <- sum(y1[-1] + y1[-n]) / 2 * delta
b <- sum(y2[-1] + y2[-n]) / 2 * delta
a / b
}
``````

Those two functions return rather similar result as we can test:

``````f1 (0.1)
# [1] -1.276167

f2 (0.1, 1000)
# [1] -1.276166
``````

Now, the only thing of interest is the limiting behaviour when `l -> 0` (in a numerical sense). Let's try

``````l <- 10 ^ (- (1:16))
# [1] 1e-01 1e-02 1e-03 1e-04 1e-05 1e-06 1e-07 1e-08 1e-09 1e-10 1e-11 1e-12
# [13] 1e-13 1e-14 1e-15 1e-16

y1 <- sapply(l, f1)
# [1] -1.2761674 -0.8698411 -0.8096179 -0.7996069 -0.7981338 -0.7979341
# [7] -0.7978877 -0.7978848 -0.7978846 -0.7978846 -0.7978846 -0.7978846
# [13] -0.7978846 -0.7978846 -0.7978846 -0.7978846

## quite a dense grid; takes some time to compute
y2 <- sapply(l, f2, n = 1e+6)
# [1] -1.2761674 -0.8698411 -0.8096179 -0.7996071 -0.7981158 -0.7979137
# [7] -0.7978877 -0.7978834 -0.7978816 -0.7978799 -0.7978783 -0.7978767
# [13] -0.7978750 -0.7978734 -0.7978717 -0.7978700
``````

Now, it looks like there is a limit toward around `-0.7978` as `l -> 0`.

Note, the `-0.8676396` you got is actually about `f1(0.01)` or `f2(0.01, 1e+6)`.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download