#
# 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<maxiter)){
iter <- iter+1
dm <- (d0+d1)/2
p0 <- prfn(d0,n,m,distrib)
pm <- prfn(dm,n,m,distrib)
p1 <- prfn(d1,n,m,distrib)
if ((p1-p)*(p-pm) > 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))
}