Nothing
saddle = function(x,lambda){
d = max(lambda)
lambda = lambda/d
x = x/d
k0 = function(zeta) -sum(log(1-2*zeta*lambda))/2
kprime0 = function(zeta) sapply(zeta, function(zz) sum(lambda/(1-2*zz*lambda)))
kpprime0 = function(zeta) 2*sum(lambda^2/(1-2*zeta*lambda)^2)
n = length(lambda)
if (any(lambda < 0)) {
lmin = max(1/(2 * lambda[lambda < 0])) * 0.99999
} else if (x>sum(lambda)){
lmin = -0.01
} else {
lmin = -length(lambda)/(2*x)
}
lmax = min(1/(2*lambda[lambda>0]))*0.99999
hatzeta = uniroot(function(zeta) kprime0(zeta) - x, lower = lmin, upper = lmax, tol = 1e-08)$root
w = sign(hatzeta)*sqrt(2*(hatzeta*x-k0(hatzeta)))
v = hatzeta*sqrt(kpprime0(hatzeta))
if(abs(hatzeta)<1e-4){
return(NA)
} else{
return( pnorm(w+log(v/w)/w, lower.tail=FALSE) )
}
}
Sadd.pval = function(Q.all,lambda){
sad = rep(1,length(Q.all))
if(sum(Q.all>0)>0){
sad[Q.all>0] = sapply(Q.all[Q.all>0],saddle,lambda=lambda)
}
id = which(is.na(sad))
if(length(id)>0){
sad[id] = Liu.pval(Q.all[id], lambda)
}
return(sad)
}
Liu.pval = function(Q, lambda){
c1 = rep(0,4); for(i in 1:4){ c1[i] = sum(lambda^i) }
muQ = c1[1]; sigmaQ = sqrt(2 *c1[2])
s1 = c1[3]/c1[2]^(3/2); s2 = c1[4]/c1[2]^2
if(s1^2 > s2){
a = 1/(s1 - sqrt(s1^2 - s2)); d = s1 *a^3 - a^2; l = a^2 - 2*d
} else {
l = 1/s2; a = sqrt(l); d = 0
}
muX = l+d; sigmaX = sqrt(2)*a
Q.Norm = (Q - muQ)/sigmaQ
Q.Norm1 = Q.Norm*sigmaX + muX
pchisq(Q.Norm1, df = l,ncp=d, lower.tail=FALSE)
}
Liu.qval.mod = function(pval, lambda){
c1 = rep(0,4)
c1[1] = sum(lambda); c1[2] = sum(lambda^2)
c1[3] = sum(lambda^3); c1[4] =sum(lambda^4)
muQ = c1[1]; sigmaQ = sqrt(2 *c1[2])
s1 = c1[3]/c1[2]^(3/2); s2 = c1[4]/c1[2]^2
beta1= sqrt(8)*s1; beta2 = 12*s2; type1 = 0
if(s1^2 > s2){
a = 1/(s1 - sqrt(s1^2 - s2)); d = s1 *a^3 - a^2; l = a^2 - 2*d
} else {
type1 = 1; l = 1/s2; a = sqrt(l); d = 0
}
muX = l+d; sigmaX = sqrt(2) *a
df = l
q.org = qchisq(pval,df=df,lower.tail=FALSE)
(q.org - df)/sqrt(2*df)*sigmaQ + muQ
}
KAT.pval <- function(Q.all, lambda, acc=1e-27,lim=1e6){
pval = rep(0, length(Q.all))
i1 = which(is.finite(Q.all))
for(i in i1){
tmp = davies(Q.all[i],lambda,acc=acc,lim=lim); pval[i] = tmp$Qq
if((tmp$ifault>0)|(pval[i]<=0)|(pval[i]>=1)) pval[i] = Sadd.pval(Q.all[i],lambda)
}
return(pval)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.