[Rd] Getting high precision values from qnorm in the tail

[Spencer Graves]

rtfm:  help('qnorm') identifies arguments you overlooked.  1-x generates 
roundoff error.  Try the following:


x <- 10^(-(1:10))
qx <- qnorm(x)
q1x <- qnorm(1-x)
qlx <- qnorm(x, lower=FALSE)

 > cbind(x, qx, q1x, qlx, qx.1x=qx+q1x, qx.lx=qx+qlx)
           x        qx      q1x      qlx         qx.1x qx.lx
  [1,] 1e-01 -1.281552 1.281552 1.281552  0.000000e+00     0
  [2,] 1e-02 -2.326348 2.326348 2.326348  0.000000e+00     0
  [3,] 1e-03 -3.090232 3.090232 3.090232  0.000000e+00     0
  [4,] 1e-04 -3.719016 3.719016 3.719016  2.842171e-14     0
  [5,] 1e-05 -4.264891 4.264891 4.264891  1.014300e-12     0
  [6,] 1e-06 -4.753424 4.753424 4.753424 -5.809575e-12     0
  [7,] 1e-07 -5.199338 5.199338 5.199338  9.784440e-11     0
  [8,] 1e-08 -5.612001 5.612001 5.612001 -8.692842e-10     0
  [9,] 1e-09 -5.997807 5.997807 5.997807  4.593952e-09     0
[10,] 1e-10 -6.361341 6.361341 6.361341 -1.270663e-08     0


hope this helps.  Spencer Graves


On 2017-04-16 6:30 AM, Sheldon Maze wrote:
> Hello All
>
> I am looking for high precision values for the normal distribution in the
> tail,(1e-10 and 1 - 1e-10) as the R package that I am using sets any number
> which is out of this range to these values and then calls the qnorm and qt
> function.
>
> What I have noticed is that the qnorm implementation in R is not symmetric
> when looking at the tails. This is quite surprising to me, as it is well
> known that this distribution is symmetric, and I have seen implementations
> in other languages that are symmetric. I have checked the qt function and
> it is also not symmetric in the tails.
>
> Here are the results from the qnorm function:
>
> x       qnorm(x)                qnorm(1-x)              qnorm(1-x) +
> qnorm(x)
> 1e-2    -2.3263478740408408     2.3263478740408408      0.0 (i.e < machine
> epsilon)
> 1e-3    -3.0902323061678132     3.0902323061678132      0.0 (i.e < machine
> epsilon)
> 1e-4    -3.71901648545568       3.7190164854557084
>   2.8421709430404007e-14
> 1e-5    -4.2648907939228256     4.2648907939238399
>   1.014299755297543e-12
> 1e-10   -6.3613409024040557     6.3613408896974208
>   -1.2706634855419452e-08
>
> It is quite clear that at a value of x close to 0 or 1, this function
> breaks down. Yes, in "normal" use this isn't a problem, but I am looking at
> fringe cases and multiplying small probabilities by very large values, in
> which case the error (1e-08) becomes a large value.
>
> Note: I have tried this with 1-x and with entering the actual number
> 0.00001 and 0.99999 and the accuracy issue is still there.
>
> The questions
>
> Firstly, is this a known problem with the qnorm and qt implementations? I
> could not find anything in the documentation, the algorithm is supposed to
> be accurate 16 digits for p values from 10^-314 as described in the
> Algorithm AS 241 paper.
>
> Quote from R doc for qnorm:
>
> "Wichura, M. J. (1988) Algorithm AS 241: The percentage points of the
> normal distribution. Applied Statistics, 37, 477–484.
>
> which provides precise results up to about 16 digits."
>
> If the R code implements the 7 digit version, why does it claim 16 digits?
> Or is it "accurate" but the original algorithm is not symmetric and wrong?
>
> If R does implement both versions of Algorithm AS 241 can I turn the 16
> digit version on?
>
> Or, is there a more accurate version of qnorm in R? Or, another solution to
> my problem where I need high precision in the tails for quantile functions.
>
> On a side note, I also have this issue with the qt distribution, at a
> similar level of precision, it is not symmetric, nor precise, but I have
> not investigated it yet. Also, I've posted this question on stack overflow:
> http://stackoverflow.com/questions/43362644/getting-high-precision-values-from-qnorm-in-the-tail
>
>
> R version:
>> version
> platform       x86_64-w64-mingw32
> arch           x86_64
> os             mingw32
> system         x86_64, mingw32
> status
> major          3
> minor          3.2
> year           2016
> month          10
> day            31
> svn rev        71607
> language       R
> version.string R version 3.3.2 (2016-10-31)
> nickname       Sincere Pumpkin Patch
>
>
> Kind regards
>
> Sheldon Maze
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel