##################################################################
# FUNCTION integral: returns value of numerical integration
##################################################################

integral <- function(f, lower, upper, ...)
{
	results <- integrate(f, lower, upper, ...)
	if(results$message != "normal termination")
		results$message
	else results$integral
}


##################################################################
# FUNCTION integralplussubdiv: returns value of numerical integration
# varying subdivisions
##################################################################

integral.plussubdiv <- function(f, lower, upper, subdiv = 100, plussubdiv = 50, maxit = 100, ...)
{
	results <- integrate(f, lower, upper, subdivisions = subdiv, ...)
	it <- 0
	while((results$message != "normal termination") && ((it <- it + 1) < maxit)) {
		subdiv <- subdiv + plussubdiv
		results <- integrate(f, lower, upper, subdivisions = subdiv, ...)
	}
	if(results$message != "normal termination")
		results$message
	else results$integral
}


##################################################################
# FUNCTION pchisq2: returns value of joint distribution function of 
# two sum of chi-squared random variables:
# Pr( chi2(dfvec[1])<=evec[1], chi2(dfvec[2])<=evec[2] )
############################################################################

pchisq2 <- function(dfvec, evec, ...)
{
	e1 <- evec[1]
	e2 <- evec[2]
	df1 <- dfvec[1]
	df2 <- dfvec[2] - dfvec[1]
	inner.int <- function(x, e2, df1, df2)
	{
		pchisq(e2 - x, df2) * dchisq(x, df1)
	}
	integral.plussubdiv(f = inner.int, lower = 0, upper = e1, e2 = e2, df1 = df1, df2 = df2)
}


##################################################################
# FUNCTION pchisq3: returns value of joint distribution function of 
# three sum of chi-squared random variables
# Pr( chi2(dfvec[1])<=evec[1], chi2(dfvec[2])<=evec[2], chi2(dfvec[3])<=evec[3] )
##################################################################

pchisq3 <- function(dfvec, evec, ...)
{
	e1 <- evec[1]
	e2 <- evec[2]
	e3 <- evec[3]
	df1 <- dfvec[1]
	df2 <- dfvec[2] - dfvec[1]
	df3 <- dfvec[3] - dfvec[2]
	inner.int <- function(x1, e2, e3, df1, df2, df3)
	{
		sapintegral <- function(x1, df1, df2, df3, e2, e3)
		{
			integral.plussubdiv(f = function(x2, x1, df2, df3, e3)
			{
				pchisq(e3 - x2 - x1, df3) * dchisq(x2, df2)
			}
			, lower = 0, upper = e2 - x1, x1 = x1, df2 = df2, df3 = df3, e3 = e3) * dchisq(x1, df1)
		}
		sapply(x1, sapintegral, df1 = df1, df2 = df2, df3 = df3, e2 = e2, e3 = e3)
	}
	integral.plussubdiv(f = inner.int, lower = 0, upper = e1, df1 = df1, df2 = df2, df3 = df3, e2 = e2, e3 = e3)
}


######### example: Table 1 ###########

# DF(2,6,10)
dfvec <- c(2,6,10) 
beta <- 0.01
divector01 <- qchisq(1-beta, df=dfvec)
list(alpha=1 - jcdfport(dfvec, divector01), beta=beta, divec=divector01)
1-pchisq3(dfvec, divector01)

# DF(12,18)
dfvec <- c(12,18) 
beta <- 0.10
divector10 <- qchisq(1-beta, df=dfvec)
list(alpha=1 - jcdfport(dfvec, divector10), beta=beta, divec=divector10)
1-pchisq2(dfvec, divector10)