[Rd] Calculation of e^{z^2/2} for a normal deviate z

[Gabriel Becker]

Hi Jing,

Peter pointed out how you can, more or less, get numbers for this, and he's
absolutely right. At the risk of giving unsolicited advice, though, Im
don't think you *should* in this case.

Someone on this list with more applied Statistics or Statistical Genetics
experience can correct me if I'm wrong, but I have always been under the
impression that "pvalues" that small  (and honestly ones that are not
nearly that small) are not meaningful as actual numbers and should never be
used as such. Essentially, AFAIK, a p-value of 1e-300 just means "really
small"/"really unlikely", the same as a pvalue of 1e-280 does (even though
the second pvalue is *20 orders of magnitude* larger). Like, how
*precicely* correct
would the distributional assumptions the test is making have to be for that
-300 to be meaningful?

Note that if you have putative pvalues in the 1e-300 range, a full order of
magnitude increase in the pvalue (to 1e-299) only corresponds to only a
reduction of 0.06 in the quantile value. Is that amount of measurement
difference likely to be meaningful in a genomics setting or is it probably
machine noise? (my impression is the later).

> log(1e-300)

[1] -690.7755

> log(1e-299)

[1] -688.4729

> qnorm(-691, log = TRUE)

[1] -37.05315

> qnorm(-688, log=TRUE)

[1] -36.97216

> qnorm(-688.4729, log=TRUE) - qnorm(-690.7755, log = TRUE)

[1] 0.06216021

 And all of the above ignores  the most important thing in practice, though
(sorry for burying the lede), which is the fact that the precision
tolerance of the computation calculating the pvalue is very likely *hundreds
of orders of magnitude *larger (ie less precise) than 1e-300 unless it uses
infinite precision arithmetic (which, I think, would be a very weird thing
to do for the reasons pointed out above).

Personally, whenever I See values on the order of 1e-K where K is bigger
than like 20 on the ouside , I just see rounding errors around 0, not
numbers that are, e.g., meaningfully comparable or sortable or summable,
etc.

So R will give you numbers for this, as Peter showed you, but personally
remain pretty skeptical that they will actually be useful for what you want
to do with them i this case.

Best,
~G

On Fri, Jun 21, 2019 at 9:32 AM peter dalgaard <pdalgd using gmail.com> wrote:

> You may want to look into using the log option to qnorm
>
> e.g., in round figures:
>
> > log(1e-300)
> [1] -690.7755
> > qnorm(-691, log=TRUE)
> [1] -37.05315
> > exp(37^2/2)
> [1] 1.881797e+297
> > exp(-37^2/2)
> [1] 5.314068e-298
>
> Notice that floating point representation cuts out at 1e+/-308 or so. If
> you want to go outside that range, you may need explicit manipulation of
> the log values. qnorm() itself seems quite happy with much smaller values:
>
> > qnorm(-5000, log=TRUE)
> [1] -99.94475
>
> -pd
>
> > On 21 Jun 2019, at 17:11 , jing hua zhao <jinghuazhao using hotmail.com>
> wrote:
> >
> > Dear Rui,
> >
> > Thanks for your quick reply -- this allows me to see the bottom of this.
> I was hoping we could have a handle of those p in genmoics such as 1e-300
> or smaller.
> >
> > Best wishes,
> >
> >
> > Jing Hua
> >
> > ________________________________
> > From: Rui Barradas <ruipbarradas using sapo.pt>
> > Sent: 21 June 2019 15:03
> > To: jing hua zhao; r-devel using r-project.org
> > Subject: Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
> >
> > Hello,
> >
> > Well, try it:
> >
> > p <- .Machine$double.eps^seq(0.5, 1, by = 0.05)
> > z <- qnorm(p/2)
> >
> > pnorm(z)
> > # [1] 7.450581e-09 1.228888e-09 2.026908e-10 3.343152e-11 5.514145e-12
> > # [6] 9.094947e-13 1.500107e-13 2.474254e-14 4.080996e-15 6.731134e-16
> > #[11] 1.110223e-16
> > p/2
> > # [1] 7.450581e-09 1.228888e-09 2.026908e-10 3.343152e-11 5.514145e-12
> > # [6] 9.094947e-13 1.500107e-13 2.474254e-14 4.080996e-15 6.731134e-16
> > #[11] 1.110223e-16
> >
> > exp(z*z/2)
> > # [1] 9.184907e+06 5.301421e+07 3.073154e+08 1.787931e+09 1.043417e+10
> > # [6] 6.105491e+10 3.580873e+11 2.104460e+12 1.239008e+13 7.306423e+13
> > #[11] 4.314798e+14
> >
> >
> > p is the smallest possible such that 1 + p != 1 and I couldn't find
> > anything to worry about.
> >
> >
> > R version 3.6.0 (2019-04-26)
> > Platform: x86_64-pc-linux-gnu (64-bit)
> > Running under: Ubuntu 19.04
> >
> > Matrix products: default
> > BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
> > LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
> >
> > locale:
> >  [1] LC_CTYPE=pt_PT.UTF-8       LC_NUMERIC=C
> >  [3] LC_TIME=pt_PT.UTF-8        LC_COLLATE=pt_PT.UTF-8
> >  [5] LC_MONETARY=pt_PT.UTF-8    LC_MESSAGES=pt_PT.UTF-8
> >  [7] LC_PAPER=pt_PT.UTF-8       LC_NAME=C
> >  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> > [11] LC_MEASUREMENT=pt_PT.UTF-8 LC_IDENTIFICATION=C
> >
> > attached base packages:
> > [1] stats     graphics  grDevices utils     datasets  methods
> > [7] base
> >
> > other attached packages:
> >
> > [many packages loaded]
> >
> >
> > Hope this helps,
> >
> > Rui Barradas
> >
> > �s 15:24 de 21/06/19, jing hua zhao escreveu:
> >> Dear R-developers,
> >>
> >> I am keen to calculate exp(z*z/2) with z=qnorm(p/2) and p is very
> small. I wonder if anyone has experience with this?
> >>
> >> Thanks very much in advance,
> >>
> >>
> >> Jing Hua
> >>
> >>       [[alternative HTML version deleted]]
> >>
> >> ______________________________________________
> >> R-devel using r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-devel
> >>
> >
> >       [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-devel using r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Office: A 4.23
> Email: pd.mes using cbs.dk  Priv: PDalgd using gmail.com
>
> ______________________________________________
> R-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

	[[alternative HTML version deleted]]