###################################################
### chunk number 1: 
###################################################
library(BGmix)
data(ybar,ss)


###################################################
### chunk number 2: 
###################################################
outdir <- BGmix(ybar, ss, nreps=c(8,8),niter=1000,nburn=1000)


###################################################
### chunk number 3: 
###################################################
params <- ccParams(file=outdir)


###################################################
### chunk number 4: 
###################################################
plotBasic(params,ybar,ss)


###################################################
### chunk number 5: 
###################################################
par(mfrow=c(1,2))
fdr <- calcFDR(params)
plotFDR(fdr)


###################################################
### chunk number 6: 
###################################################
pred <- ccPred(file=outdir)


###################################################
### chunk number 7: 
###################################################
par(mfrow=c(1,2))
hist(pred$pval.ss.mix[,1])
hist(pred$pval.ss.mix[,2])


###################################################
### chunk number 8: 
###################################################
par(mfrow=c(2,3))
plotPredChecks(pred$pval.ybar.mix2,params$pc,probz=0.8)


###################################################
### chunk number 9: 
###################################################
nreps <- c(8,8)
outdir2 <- BGmix(ybar, ss, nreps=nreps, jstar=-1, niter=1000,nburn=1000)
params2 <- ccParams(outdir2)
res2 <-  ccTrace(outdir2)


###################################################
### chunk number 10: 
###################################################
tpp.res <- TailPP(res2, nreps=nreps, params2, p.cut = 0.7, plots  = F)


###################################################
### chunk number 11: 
###################################################
FDR.res <- FDRforTailPP(tpp.res$tpp, a1 = params2$maa[1], a2
= params2$maa[2], n.rep1=nreps[1], n.rep2=nreps[2], p.cut = 0.8)
pi0 <- EstimatePi0(tpp.res$tpp, tpp.res$pp0)


###################################################
### chunk number 12: 
###################################################
par(mfrow=c(1,2))
histTailPP(tpp.res)
FDRplotTailPP(tpp.res)


###################################################
### chunk number 13: 
###################################################
unlink("run.1", recursive=TRUE)
unlink("run.2", recursive=TRUE)