vrsn <- "20120804"
Dirname <- "Inference"
FileName <- "RCTvalidateInference"

RCTvalidateOutputDirname <- paste(RCTdirnameV,Dirname,sep="")
if(!file.exists(RCTvalidateOutputDirname))
	dir.create(RCTvalidateOutputDirname)
RCTvalidateOutputDirname <- paste(RCTvalidateOutputDirname,"/",sep="")
cat("#####", FileName, "-", vrsn, "\n")
if (!exists("runRCTvalidateTests")) {
	runRCTvalidateTests <- 1:38
} else {
	cat("   ## Tests ", runRCTvalidateTests, "\n")
}

gstMLEle <- function (M, X, dsn, m, x, theta) {
	bnd <- seqBoundary (dsn, "X")
	Nj <- seqExtract(dsn,"sample.size")
	J <- length(Nj)
	if (J==1) antmOK <- TRUE
	else antmOK <- all(attr(bnd,"no.stopping")[-J])	
	if (length(M) != length(X)) {
		if (length(M)==1) M <- rep(M,length(S))
		else if (length(S)==1) S <- rep(S,length(M))
		else stop("M, S should have same length")
	} 
	if (!all(c(m,M) %in% 1:J)) stop("m, M must agree with dsn")
	if ((x > bnd[m,1] & x < bnd[m,2]) | (x > bnd[m,3] & x < bnd[m,4])) stop("(m,x) not a valid sequential observation")
	if (any((X > bnd[M,1] & X < bnd[M,2]) | (X > bnd[M,3] & X < bnd[M,4])))
		stop("not all (M,X) are valid sequential observations")
	if (antmOK) antm <- ((M==m) & (X <= x)) | ((M < m) & (X < bnd[M,1])) | ((M > m) & (x > bnd[m,4]))
	else antm <- rep(NA,length(M))
	mle <- X <= x
	lr <- sqrt(Nj[M])*(X - theta) <= sqrt(Nj[m])*(x - theta)
	list(meanOrder=mle,timeOrder=antm,LROrder=lr)
}

seqCheckQuantilesBAM <- function (dsn, theta=NULL, probs=NULL, Nsimul=10000, ...) {
	if (!is.null(theta)) zOC <- seqOC(dsn, theta=theta[1])
	else zOC <- seqOC(dsn, power=seqExtract(dsn, "power"))
	theta <- zOC$theta
	zQ <- qSeq(dsn, prob=probs, theta=theta)
	antmOK <- all(!is.na(zQ$analysis.index))
	x <- rSeq(Nsimul, dsn, theta=theta, ...)
	expMLE <- c(theta, zOC$expTheta, mean(x$thetaHat), sqrt(var(x$thetaHat)))
	expMLE <- c(expMLE, sqrt(Nsimul) * (expMLE[3] - expMLE[2]) / expMLE[4])
	expMLE <- c(expMLE, 2 * pnorm(- abs(expMLE[5])))
	names(expMLE) <- c("Theta", "expTheta", "meanThetaHat", "sdThetaHat", "Z", "Pval")
	simQmean <- simQLR <- simQtime <- cbind(Prob=probs, EmpProb=0, Z=0, Pval=0)
	for (j in 1:length(probs)) {
		simQmean[j,2] <- mean(gstMLEle (x$An.Time, x$thetaHat, dsn, J, zQ$qObsMean[j], theta)$meanOrder)
		if (antmOK) simQtime[j,2] <- mean(gstMLEle (x$An.Time, x$thetaHat, dsn, zQ$analysis.index[j], zQ$observed[j], theta)$timeOrder)
		else simQtime[j,2] <- NA
		simQLR[j,2] <- mean(gstMLEle (x$An.Time, x$thetaHat, dsn, J, zQ$qLR[j], theta)$LROrder)
	}
	simQmean[,3] <- (simQmean[,2] - probs) / sqrt(probs * (1 - probs) / Nsimul)
	simQmean[,4] <- 2 * pnorm(-abs(simQmean[,3]))
	simQtime[,3] <- (simQtime[,2] - probs) / sqrt(probs * (1 - probs) / Nsimul)
	simQtime[,4] <- 2 * pnorm(-abs(simQtime[,3]))
	simQLR[,3] <- (simQLR[,2] - probs) / sqrt(probs * (1 - probs) / Nsimul)
	simQLR[,4] <- 2 * pnorm(-abs(simQLR[,3]))
	rr <- rep(probs,each=length(probs))
	cc <- rep(probs,length(probs))
	V <- solve((ifelse(rr < cc, rr, cc) - probs %o% probs) / Nsimul)
	statQ <- c(((simQmean[,2] - probs) %*% V) %*% (simQmean[,2] - probs),
		((simQLR[,2] - probs) %*% V) %*% (simQLR[,2] - probs),
		((simQtime[,2] - probs) %*% V) %*% (simQtime[,2] - probs))
	statQ <- cbind(statQ, 1-pchisq(as.vector(statQ),length(probs)))
	dimnames(statQ) <- list(c("meanQ","LRQ","timeQ"),c("Chi2stat","Pval"))
	list(expMLE=expMLE, statQ=statQ, simQmean=simQmean, simQtime=simQtime, simQLR=simQLR, design=dsn, Nsimul=Nsimul)
}

seqCheckConsistentInference <- function (dsn, j, thetaHat, dsnInfrnc=NULL) {
	if (is.null(dsnInfrnc)) dsnInfrnc <- seqInference(dsn, analysis.index=j[1], obs=thetaHat[1], conf.level= 0.95)
	else if(is.seqInference(dsnInfrnc)) dsnInfrnc <- dsnInfrnc[1,]
	else stop("dsnInfrnc needs to be a seqInference object")
	c(BAM = thetaHat - seqOC(dsn, theta= dsnInfrnc$BAM)$expTheta,
		MUE.m = 0.5 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$MUE.meanOrder)$cdf.meanOrder,
		MUE.l = 0.5 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$MUE.LROrder)$cdf.LROrder,
		MUE.t = 0.5 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$MUE.timeOrder)$cdf.timeOrder,
		CIlo.m = 0.975 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$CIlo.meanOrder)$cdf.meanOrder,
		CIhi.m = 0.025 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$CIhi.meanOrder)$cdf.meanOrder,
		CIlo.l = 0.975 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$CIlo.LROrder)$cdf.LROrder,
		CIhi.l = 0.025 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$CIhi.LROrder)$cdf.LROrder,
		CIlo.t = 0.975 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$CIlo.timeOrder)$cdf.timeOrder,
		CIhi.t = 0.025 - pSeq(dsn, j, thetaHat, theta= dsnInfrnc$CIhi.timeOrder)$cdf.timeOrder
		)
}

seqSimCoverageBias <- function (dsn, theta, Nsimul= 2000, conf.level=0.95, ...) {
	theta <- theta[1]
	x <- rSeq(Nsimul, dsn, theta=theta, ...)
	zz <- seqInference(z, analysis.index=x$An.Time, observed=x$thetaHat, conf.level=conf.level)
	cvrg <- c(theta, mean(zz$CIlo.meanOrder < theta & zz$CIhi.meanOrder > theta),
		mean(zz$CIlo.LROrder < theta & zz$CIhi.LROrder > theta),
		mean(zz$CIlo.timeOrder < theta & zz$CIhi.timeOrder > theta))
	avglng <- c(theta, mean(zz$CIlo.meanOrder - zz$CIhi.meanOrder),
		mean(zz$CIlo.LROrder - zz$CIhi.LROrder),
		mean(zz$CIlo.timeOrder - zz$CIhi.timeOrder))
	bias <- c(theta, mean(zz$MLE), mean(zz$BAM), mean(zz$RBadj), mean(zz$MUE.meanOrder), 
		mean(zz$MUE.LROrder), mean(zz$MUE.timeOrder))
	bias[-1] <- bias[-1] - theta
	rmse <- c(theta, mean((zz$MLE-theta)^2), mean((zz$BAM-theta)^2), mean((zz$RBadj-theta)^2), 
		mean((zz$MUE.meanOrder-theta)^2), mean((zz$MUE.LROrder-theta)^2), mean((zz$MUE.timeOrder-theta)^2))
	rmse[-1] <- sqrt(rmse[-1])
	names(cvrg) <- names(avglng) <- c("Theta", "meanOrder", "LROrder", "timeOrder")
	PvalCvrg <- cvrg
	PvalCvrg[2:4] <- 2 * pnorm(-abs((PvalCvrg[2:4] - conf.level) / sqrt(conf.level * (1 - conf.level) / Nsimul)))
	names(bias) <- names(rmse) <- c("Theta", "MLE", "BAM", "RBadj", "MUE.meanOrder", "MUE.LROrder", "MUE.timeOrder")
	list(Coverage=cvrg, PvalCvrg=PvalCvrg, AvgLength=avglng, Bias=bias, RMSE=rmse, Nsimul=Nsimul, design=dsn)
}

####################
TestName <- "One Sample Normal "
####################

rslt <- NULL
testLbls <- (1:17)[(1:17) %in% runRCTvalidateTests]

# Tests 1.1-1.4
test <- 1
if (test %in% runRCTvalidateTests) {
	testLbls <- c(testLbls,paste(test,1:4,sep="."))
	dsnName <- "OBF1s.025.5"
	J <- 5
	P <- 1
	alpha <- 0.025
	beta <- 0.975
	z <- seqDesign (arms=1,sample.size=100,nbr.analyses=J,alpha=alpha,beta=beta,P=P,suppressErrors=T,EPSILON=1e-8)
	set.seed(test)
	theta <- seqOC(z, power=runif(1))$theta
	check <- TRUE
	x <- rSeq(5, z, theta=theta, var.true=T, sufficient=T)
	zz <- seqInference(z, analysis.index=x$An.Time, observed=x$thetaHat)
	for (j in 1:5)
		check <- check && 
			all(abs(seqCheckConsistentInference (z, zz$analysis.index[j], zz$observed[j], zz[j,])) < 1e-8, na.rm=TRUE)
	###### test 1.1 consistency of inference with operating characteristics ######
	rslt <- c(rslt,check)
	zz1 <- seqCheckQuantilesBAM (z, theta=theta, probs=c(1:4,5*(1:19),96:99)/100, Nsimul=10000, 
		var.true=T, sufficient=T)
	###### test 1.2 accuracy of quantiles                                   ######
	rslt <- c(rslt,all(zz1$statQ[,2] > .05, na.rm=TRUE))
	###### test 1.3 accuracy of expected MLE                                ######
	rslt <- c(rslt,zz1$expMLE["Pval"] > .05)
	zz2 <- seqSimCoverageBias (z, theta=theta, var.true=T, sufficient=T)
	###### test 1.4 coverage accuracy of confidence intervals               ######
	rslt <- c(rslt,all(zz2$PvalCvrg[2:4] > .05, na.rm=TRUE))
}

# Tests 2.1-2.4
test <- 2
if (test %in% runRCTvalidateTests) {
	testLbls <- c(testLbls,paste(test,1:4,sep="."))
	dsnName <- "Poc1s.025.3"
	J <- 3
	P <- 0.5
	alpha <- 0.025
	beta <- 0.975
	z <- seqDesign (arms=1,sample.size=100,nbr.analyses=J,alpha=alpha,beta=beta,P=P,suppressErrors=T,EPSILON=1e-8)
	set.seed(test)
	theta <- seqOC(z, power=runif(1))$theta
	check <- TRUE
	x <- rSeq(5, z, theta=theta, var.true=T, sufficient=T)
	zz <- seqInference(z, analysis.index=x$An.Time, observed=x$thetaHat)
	for (j in 1:5)
		check <- check && 
			all(abs(seqCheckConsistentInference (z, zz$analysis.index[j], zz$observed[j], zz[j,])) < 1e-8, na.rm=TRUE)
	###### test 2.1 consistency of inference with operating characteristics ######
	rslt <- c(rslt,check)
	zz1 <- seqCheckQuantilesBAM (z, theta=theta, probs=c(1:4,5*(1:19),96:99)/100, Nsimul=10000, 
		var.true=T, sufficient=T)
	###### test 2.2 accuracy of quantiles                                   ######
	rslt <- c(rslt,all(zz1$statQ[,2] > .05, na.rm=TRUE))
	###### test 2.3 accuracy of expected MLE                                ######
	rslt <- c(rslt,zz1$expMLE["Pval"] > .05)
	zz2 <- seqSimCoverageBias (z, theta=theta, var.true=T, sufficient=T)
	###### test 2.4 coverage accuracy of confidence intervals               ######
	rslt <- c(rslt,all(zz2$PvalCvrg[2:4] > .05, na.rm=TRUE))
}

# Tests 3.1-3.4
test <- 3
if (test %in% runRCTvalidateTests) {
	testLbls <- c(testLbls,paste(test,1:4,sep="."))
	dsnName <- "PT1s.025.5"
	J <- 5
	P <- 0.5
	alpha <- c(0.2,0.025)
	beta <- c(0.975,0.8)
	z <- seqDesign (arms=1,sample.size=100,nbr.analyses=J,alpha=alpha,beta=beta,P=P,suppressErrors=T,EPSILON=1e-8)
	set.seed(test)
	theta <- seqOC(z, power=runif(1))$theta
	check <- TRUE
	x <- rSeq(5, z, theta=theta, var.true=T, sufficient=T)
	zz <- seqInference(z, analysis.index=x$An.Time, observed=x$thetaHat)
	for (j in 1:5)
		check <- check && 
			all(abs(seqCheckConsistentInference (z, zz$analysis.index[j], zz$observed[j], zz[j,])) < 1e-8, na.rm=TRUE)
	###### test 3.1 consistency of inference with operating characteristics ######
	rslt <- c(rslt,check)
	zz1 <- seqCheckQuantilesBAM (z, theta=theta, probs=c(1:4,5*(1:19),96:99)/100, Nsimul=10000, 
		var.true=T, sufficient=T)
	###### test 3.2 accuracy of quantiles                                   ######
	rslt <- c(rslt,all(zz1$statQ[,2] > .05, na.rm=TRUE))
	###### test 3.3 accuracy of expected MLE                                ######
	rslt <- c(rslt,zz1$expMLE["Pval"] > .05)
	zz2 <- seqSimCoverageBias (z, theta=theta, var.true=T, sufficient=T)
	###### test 3.4 coverage accuracy of confidence intervals               ######
	rslt <- c(rslt,all(zz2$PvalCvrg[2:4] > .05, na.rm=TRUE))
}

# Tests 4.1-4.4
test <- 4
if (test %in% runRCTvalidateTests) {
	testLbls <- c(testLbls,paste(test,1:4,sep="."))
	dsnName <- "UF2s.5"
	J <- 5
	P <- c(0.5,1)
	alpha <- c(0.2,0.025)
	beta <- c(0.975,0.8)
	z <- seqDesign (arms=1,sample.size=100,nbr.analyses=J,alpha=alpha,beta=beta,P=P,suppressErrors=T,
		test.type="two.sided",early.stopping="both",EPSILON=1e-8)
	set.seed(test)
	theta <- seqOC(z, power=runif(1))$theta
	check <- TRUE
	x <- rSeq(5, z, theta=theta, var.true=T, sufficient=T)
	zz <- seqInference(z, analysis.index=x$An.Time, observed=x$thetaHat)
	for (j in 1:5)
		check <- check && 
			all(abs(seqCheckConsistentInference (z, zz$analysis.index[j], zz$observed[j], zz[j,])) < 1e-8, na.rm=TRUE)
	###### test 4.1 consistency of inference with operating characteristics ######
	rslt <- c(rslt,check)
	zz1 <- seqCheckQuantilesBAM (z, theta=theta, probs=c(1:4,5*(1:19),96:99)/100, Nsimul=10000, 
		var.true=T, sufficient=T)
	###### test 4.2 accuracy of quantiles                                   ######
	rslt <- c(rslt,all(zz1$statQ[,2] > .05, na.rm=TRUE))
	###### test 4.3 accuracy of expected MLE                                ######
	rslt <- c(rslt,zz1$expMLE["Pval"] > .05)
	zz2 <- seqSimCoverageBias (z, theta=theta, var.true=T, sufficient=T)
	###### test 4.4 coverage accuracy of confidence intervals               ######
	rslt <- c(rslt,all(zz2$PvalCvrg[2:4] > .05, na.rm=TRUE))
}

if (is.null(rslt)) {cat("     ",TestName,": No tests performed\n")
} else if (all(rslt)) {cat("     ",TestName,": All",length(rslt),"tests  PASS\n")
} else cat("     ",TestName,": failing test numbers ", testLbls[!rslt],"\n")


rm(FileName,TestName,Dirname,RCTvalidateOutputDirname,rslt,check,z,zz,zz1,zz2,dsnName,J,P,alpha,beta,test,testLbls,x)