# # Functions for HDRs and for Behrens' distribution # hdrsize <- function(p,n,m,distrib){ # Uses algorithm described in Jackson (1974) if (p > 1) p <- p/100 precision <- 0.0001 maxiter <- 100 pp <- 1 d0 <- dinitial(p,n,m,distrib) d1 <- d0 iter <- 0 while (prfn(d1,n,m,distrib) > p) d1 <- d0/2 while (prfn(d0,n,m,distrib) < p) d0 <- 2*d0 while ((abs(pp)>precision)&&(iter 0) d0 <- dm if ((pm-p)*(p-p0) > 0) d1 <- dm pp <- abs(p1-p0) } return(d0) } theta0fn <- function(d,n,m,distrib){ if (distrib=="gamma") return(d/tanh(d/((n-1)))) if (distrib=="invgamma") return(d/tanh(d/((n+1)))) if (distrib=="gammafromlogs") return(d/tanh(d/n)) if (distrib=="beta") return((d^((n-1)/(m-1))-1)/(d^((n-1)/(m-1)+1)-1)) if (distrib=="f") return((m/n)*(d-d^((m+2)/(m+n)))/(d^((m+2)/(m+n))-1)) if (distrib=="ffromlogs") return((m/n)*(d-d^(m/(m+n)))/(d^(m/(m+n))-1)) } prfn <- function(d,n,m,distrib){ if (distrib=="gamma") return(pgamma(upperfn(d,n,m,distrib),n,1)- pgamma(lowerfn(d,n,m,distrib),n,1)) if (distrib=="invgamma") return(pgamma(1/lowerfn(d,n,m,distrib),n,1)- pgamma(1/upperfn(d,n,m,distrib),n,1)) if (distrib=="gammafromlogs") return(pgamma(upperfn(d,n,m,distrib),n,1)- pgamma(lowerfn(d,n,m,distrib),n,1)) if (distrib=="beta") return(pbeta(upperfn(d,n,m,distrib),n,m)- pbeta(lowerfn(d,n,m,distrib),n,m)) if (distrib=="f") return(pf(upperfn(d,n,m,distrib),n,m)- pf(lowerfn(d,n,m,distrib),n,m)) if (distrib=="ffromlogs") return(pf(upperfn(d,n,m,distrib),n,m)- pf(lowerfn(d,n,m,distrib),n,m)) } dinitial <- function(p,n,m,distrib){ if (distrib=="gamma") return((qgamma(1-p/2,n,1)-qgamma(p/2,n,1))/2) if (distrib=="invgamma") return((1/qgamma(p/2,n,1)-1/qgamma(1-p/2,n,1))/2) if (distrib=="gammafromlogs") return((qgamma(1-p/2,n,1)-qgamma(p/2,n,1))/2) if (distrib=="beta") return(qbeta(1-p/2,n,m)/qbeta(p/2,n,m)) if (distrib=="f") return(qf(1-p/2,n,m)/qf(p/2,n,m)) if (distrib=="ffromlogs") return(qf(1-p/2,n,m)/qf(p/2,n,m)) } lowerfn <- function(d0,n,m,distrib){ if (distrib=="gamma") return(theta0fn(d0,n,1,distrib)-d0) if (distrib=="invgamma") return(1/(theta0fn(d0,n,1,distrib)+d0)) if (distrib=="gammafromlogs") return(theta0fn(d0,n,1,distrib)-d0) if (distrib=="beta") return(theta0fn(d0,n,m,distrib)) if (distrib=="f") return(d0^(-1)*theta0fn(d0,n,m,distrib)) if (distrib=="ffromlogs") return(d0^(-1)*theta0fn(d0,n,m,distrib)) } upperfn <- function(d0,n,m,distrib){ if (distrib=="gamma") return(theta0fn(d0,n,1,distrib)+d0) if (distrib=="invgamma") return(1/(theta0fn(d0,n,1,distrib)-d0)) if (distrib=="gammafromlogs") return(theta0fn(d0,n,1,distrib)+d0) if (distrib=="beta") return(d0*theta0fn(d0,n,m,distrib)) if (distrib=="f") return(theta0fn(d0,n,m,distrib)) if (distrib=="ffromlogs") return(theta0fn(d0,n,m,distrib)) } hnorm <- function(p,mean=0,sd=1) return(c(qnorm((1-p)/2,mean,sd),qnorm((1+p)/2,mean,sd))) ht <- function(p,df,ncp=0) return(c(qt((1-p)/2,df,ncp),qt((1+p)/2,df,ncp))) hcauchy <- function(p,location=0,scale=1){ y <- qcauchy((1-p)/2,location,scale,log=FALSE) z <- qcauchy((1+p)/2,location,scale,log=FALSE) return(c(y,z)) } hf <- function(p,df1,df2,log=TRUE){ if (log){ d0 <- hdrsize(p,df1,df2,"ffromlogs") return(c(lowerfn(d0,df1,df2,"ffromlogs"), upperfn(d0,df1,df2,"ffromlogs"))) } else{ d0 <- hdrsize(p,df1,df2,"f") return(c(lowerfn(d0,df1,df2,"f"), upperfn(d0,df1,df2,"f"))) } } hgamma <- function(p,shape,rate=1,scale=1/rate,log=TRUE,inverse=FALSE){ if (log){ if (inverse){ d0 <- hdrsize(p,shape,rate,"gammafromlogs") return(c(1/upperfn(d0,shape,rate,"gammafromlogs"), 1/lowerfn(d0,shape,rate,"gammafromlogs"))/(4*scale)) } else{ d0 <- hdrsize(p,shape,rate,"gammafromlogs") return(4*scale*c(lowerfn(d0,shape,rate,"gammafromlogs"), upperfn(d0,shape,rate,"gammafromlogs"))) } } else{ if (inverse){ d0 <- hdrsize(p,shape,rate,"invgamma") return(c(lowerfn(d0,shape,rate,"invgamma"), upperfn(d0,shape,rate,"invgamma"))/(4*scale)) } else{ } d0 <- hdrsize(p,shape,rate,"gamma") return(4*scale*c(lowerfn(d0,shape,rate,"gamma"), upperfn(d0,shape,rate,"gamma"))) } } hchisq <- function(p,df,log=TRUE,inverse=FALSE) return(hgamma(p,df/2,2,log=log,inverse=inverse)) hbeta <- function(p,shape1,shape2,log=TRUE){ if (log){ return(shape1*hf(p,2*shape1,2*shape2)/ (shape2+shape1*hf(p,2*shape1,2*shape2))) } else{ d0 <- hdrsize(p,shape1,shape2,"beta") return(c(lowerfn(d0,shape1,shape2,"beta"), upperfn(d0,shape1,shape2,"beta"))) } } # # Behrens' distribution # gbehrens <- function(x,z,n,m,phi,degrees=TRUE){ if (degrees) phi <- pi*phi/180 k <- 1/(beta(n/2,1/2)*beta(m/2,1/2)*sqrt(n*m)) return(k*(1+(z*cos(phi)-x*sin(phi))^2/n)^(-(n+1)/2)* (1+(z*sin(phi)+x*cos(phi))^2/m)^(-(m+1)/2)) } dbehrens <- function(z,n,m,phi,degrees=TRUE){ g <- function(x) gbehrens(x,z,n,m,phi,degrees) return(integrate(g,-Inf,Inf)$value) } pbehrens <- function(z,n,m,phi,degrees=TRUE){ # Uses the function adapt from the package adapt # Parameter Large needed as adapt will not accept Inf library(adapt) Large <- 100 f <- function(x) gbehrens(x[1],x[2],n,m,phi,degrees) if (z==0) y <- 0.5 if (z > 0) y <- 0.5+adapt(ndim=2, lower=c(0,-Large),upper=c(z,Large),functn=f)$value if (z < 0) y <- 0.5-adapt(ndim=2, lower=c(0,-Large),upper=c(-z,Large),functn=f)$value if (y < 0) y <- 0 if (y > 1) y <- 1 return(y) } qbehrens <- function(p,n,m,phi,degrees=TRUE){ precision <- 0.01 x <- qnorm(p) while (pbehrens(x,n,m,phi,degrees) < p) x <- x + precision p0 <- pbehrens(x-precision,n,m,phi,degrees) p1 <- pbehrens(x,n,m,phi,degrees) return(x - precision*(p1-p)/(p1-p0)) } rbehrens <- function(k,n,m,phi,degrees=TRUE){ y <- runif(k) for (i in 1:k) y[i] <- qbehrens(y[i],n,m,phi) return(y) } hbehrens <- function(p,n,m,phi,degrees=TRUE){ y <- qbehrens((1-p)/2,n,m,phi,degrees) z <- qbehrens((1+p)/2,n,m,phi,degrees) x <- (z-y)/2 return(c(-x,x)) }