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.
______________________________________________
R-help@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.
