R help - Nov 2011 - How to increase precision to handle very low P-values

Hello, Rlisters

I have to compute p-values that are on the tail of the distribution,
P-values < 10^-20.

However, my current implementations enable one to estimate P-values up to
10^-12, or so.

A typical example is found below, where t is my critical value.

########### example - code adapted from Rassoc #######################

rho01 = 0.5
rho105 = 0.5
rho005  = 0.5
t = 8
z = 2

w0=(rho005-rho01*rho105)/(1-rho01^2)
w1=(rho105-rho01*rho005)/(1-rho01^2)
  
 
 fun1=function(t,z){
 return(pnorm((t-rho01*z)/sqrt(1-rho01^2))*dnorm(z))
 }
 fun2=function(t,z){
 return(pnorm(((t-w0*z)/w1-rho01*z)/sqrt(1-rho01^2))*dnorm(z))
 }
 fun3=function(t,z){
 return(pnorm((-t-rho01*z)/sqrt(1-rho01^2))*dnorm(z))
 }
 
asy=function(t){
z1=2*integrate(function(z){fun1(t,z)},lower=0,upper=t*(1-w1)/w0,subdivisions=1000)$value
z2=2*integrate(function(z){fun2(t,z)},lower=t*(1-w1)/w0,upper=t,subdivisions=1000)$value
z3=-2*integrate(function(z){fun3(t,z)},lower=0,upper=t,subdivisions=500)$value
return(z1+z2+z3)
}
    
pvalue <-  1-asy(t)
pvalue
 ###########################################

Using this script, my critical values can achieve values up to 7.5, or so.
Above this cutoff my P-values show up as negative values. Why's that?


Grateful for any tips.

All the best,

Alonso


--
View this message in context:
http://r.789695.n4.nabble.com/How-to-increase-precision-to-handle-very-low-P-values-tp4100250p4100250.html
Sent from the R help mailing list archive at Nabble.com.

On 23/11/2011 11:16 AM, alonso_canada wrote:
> Hello, Rlisters
>
> I have to compute p-values that are on the tail of the distribution,
> P-values<  10^-20.
>
> However, my current implementations enable one to estimate P-values up to
> 10^-12, or so.
>
> A typical example is found below, where t is my critical value.
The p* functions (pnorm, etc.) generally have a "log" argument, so
they
can return logs of probabilities.  Together with the "lower" argument,
they have a huge range.

To use them with integrate(), rescale them.  For example, to find the 
integral of pnorm() from z=-32 to -30 you could do:

maxval <- pnorm(-30, log=TRUE)
scaled <- function(x) exp( pnorm(x, log=TRUE) - maxval )
scaledintegral <- integrate(scaled, lower=-32, upper=-30)
# Log of the result:
log(scaledintegral$value) + maxval

(The answer is -457.7, so it is actually representable:  1.631957e-199.)

Duncan Murdoch

> ########### example - code adapted from Rassoc #######################
>
> rho01 = 0.5
> rho105 = 0.5
> rho005  = 0.5
> t = 8
> z = 2
>
> w0=(rho005-rho01*rho105)/(1-rho01^2)
> w1=(rho105-rho01*rho005)/(1-rho01^2)
>
>
>   fun1=function(t,z){
>   return(pnorm((t-rho01*z)/sqrt(1-rho01^2))*dnorm(z))
>   }
>   fun2=function(t,z){
>   return(pnorm(((t-w0*z)/w1-rho01*z)/sqrt(1-rho01^2))*dnorm(z))
>   }
>   fun3=function(t,z){
>   return(pnorm((-t-rho01*z)/sqrt(1-rho01^2))*dnorm(z))
>   }
>
> asy=function(t){
>
z1=2*integrate(function(z){fun1(t,z)},lower=0,upper=t*(1-w1)/w0,subdivisions=1000)$value
>
z2=2*integrate(function(z){fun2(t,z)},lower=t*(1-w1)/w0,upper=t,subdivisions=1000)$value
>
z3=-2*integrate(function(z){fun3(t,z)},lower=0,upper=t,subdivisions=500)$value
> return(z1+z2+z3)
> }
>
> pvalue<-  1-asy(t)
> pvalue
>   ###########################################
>
> Using this script, my critical values can achieve values up to 7.5, or so.
> Above this cutoff my P-values show up as negative values. Why's that?
>
>
> Grateful for any tips.
>
> All the best,
>
> Alonso
>
>
> --
> View this message in context:
http://r.789695.n4.nabble.com/How-to-increase-precision-to-handle-very-low-P-values-tp4100250p4100250.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

Interesting observation, Duncan! I investigate its effects in my case. Thanks
again for your time. all the best, Alonso

--
View this message in context:
http://r.789695.n4.nabble.com/How-to-increase-precision-to-handle-very-low-P-values-tp4100250p4104503.html
Sent from the R help mailing list archive at Nabble.com.