library(mvtnorm)
ftn_ZCLRT = function(N2,N4,N5){

	n20 = N2[1]
	n21 = N2[2]

	n40 = N4[1]
	n41 = N4[2]
	n42 = N4[3]

	n51 = N5[1]
	n52 = N5[2]

	n2 = sum(N2)
	n4 = sum(N4)
	n5 = sum(N5)

	l_likelihood = function(x1,x2){
	(n21+n41+n51)*log(x1) + (n42+n52)*log(x2) - n2*log(1+x1) - n4*log(1+2*x1+x2) - n5*log(x1+x2)
	}


	x1 = seq(1,5,0.005)
	L = length(x1)
	value = numeric(L*(L+1)/2)

	ct = 0
	for(i in 1:L){
		x2 = seq(x1[i],5,0.005)
		for(j in 1:length(x2)){
		ct = ct + 1
		value[ct] = l_likelihood(x1[i],x2[j])
		}}

	ZCLRT = 2*(max(value) - l_likelihood(1,1))
	return(ZCLRT)
}

ftn_pCLRT = function(N_NULL,ZCLRT,corr_ZRD){


	z_null = numeric(N_NULL)
	for(i in 1:N_NULL){
	u = runif(1,0,1)
	p1 = acos(corr_ZRD)/(2*pi)
	if(u <p1){
	z_null[i] = rchisq(1,2)
	} else if(u < (p1+0.5)){
	z_null[i] = rchisq(1,1)
	} else{
	z_null[i] = 0
	}
	} # end i

	out = length(z_null[z_null >  ZCLRT])/N_NULL
	return(out)
} # end function
ftn_pMAX = function(N_NULL,ZMAX,corr_ZRA,corr_ZRD,corr_ZAD){

	corr_ZAR = corr_ZRA
	corr_ZDR = corr_ZRD
	corr_ZDA = corr_ZAD

	cor_mat = matrix(c(1,corr_ZRA,corr_ZRD,corr_ZAR,1,corr_ZAD,corr_ZDR,corr_ZDA,1),nrow=3,byrow=T)
	z_null = rmvnorm(N_NULL,c(0,0,0),cor_mat)
	z_null = apply(z_null,1,function(z){ 
	max((z)) })

	out = length(z_null[z_null > ZMAX])/N_NULL
	return(out)
}