###################################################
### chunk number 1: 
###################################################
library(bayesclust, version="2.1")


###################################################
### chunk number 2: 
###################################################
data.ex1 <- matrix(rnorm(40), nrow=20)
plot(data.ex1, main="Example Data 1", pch=20)


###################################################
### chunk number 3: 
###################################################
test.ex1 <- cluster.test(data.ex1, nsim=2500, p=2, aR=0.6, k=2)
plot(test.ex1, main="Convergence Plot for test.ex1")


###################################################
### chunk number 4: 
###################################################
summary(test.ex1)


###################################################
### chunk number 5: 
###################################################
str(test.ex1)


###################################################
### chunk number 6: 
###################################################
test.ex1$data1[length(test.ex1$data1)]


###################################################
### chunk number 7: 
###################################################
p <- 2
RW <- 0.6
NSIM <- 2500
K <- 2
fname.pdf <- paste(paste("convergence_rep1to4",RW,NSIM,K,sep="_"), ".pdf", sep="")
fname.csv <- paste(paste("finalprobs",RW,NSIM,K,sep="_"), ".csv", sep="")
pdf(file=fname.pdf)
par(mfrow=c(2,2))
output.probs <- vector("numeric", 4)
for (m in 1:4) {
tmp <- cluster.test(data.ex1, nsim=NSIM, k=K, aR=RW, p=p)
output.probs[m] <- tmp$data1[length(tmp$data1)]
plot(tmp)
}
dev.off()
write.csv(cbind(REP=1:4,NSIM, RW, K, output.probs), fname.csv, quote=FALSE, row.names=FALSE)
mean.output.probs <- mean(output.probs)


###################################################
### chunk number 8: 
###################################################
test.nulld <- nulldensity(nsim=1000, n=20, k=2, p=2)
hist(test.nulld, main="Histogram of null distribution")


###################################################
### chunk number 9: 
###################################################
quantile(test.nulld$gen, 0.05)


###################################################
### chunk number 10: 
###################################################
ex1.pval <- emp2pval(test.ex1, test.nulld)
ex1.pval$pvals


###################################################
### chunk number 11: 
###################################################
pvalue <- findInterval(mean.output.probs, test.nulld$gen) + 1
pvalue <- pvalue/length(test.nulld$gen.values)
pvalue <- min(1, pvalue)
pvalue


###################################################
### chunk number 12: 
###################################################
data(cutoffs)
head(cutoffs)


###################################################
### chunk number 13: 
###################################################
data.ex2 <- matrix(rnorm(40), nrow=20)
test.ex2 <- cluster.test(data.ex2, nsim=2500, p=2, aR=0.6, k=2)
data.ex3 <- matrix(rnorm(40), nrow=20)
test.ex3 <- cluster.test(data.ex3, nsim=2500, p=2, aR=0.6, k=2)
data.ex4 <- matrix(rnorm(40), nrow=20)
test.ex4 <- cluster.test(data.ex4, nsim=2500, p=2, aR=0.6, k=2)


###################################################
### chunk number 14: 
###################################################
ex2.pval <- emp2pval(test.ex2, test.nulld)
ex3.pval <- emp2pval(test.ex3, test.nulld)
ex4.pval <- emp2pval(test.ex4, test.nulld)
fdr.test(ls(pattern="pval$"), 0.30)


###################################################
### chunk number 15: 
###################################################
pval.string <- runif(10)
fdr.seq <- 1/10 * 0.3
min(sort(pval.string) >= fdr.seq)


###################################################
### chunk number 16: 
###################################################
opt.ex1 <- cluster.optimal(data.ex1, nsim=2500)
plot(opt.ex1)


###################################################
### chunk number 17: 
###################################################
opt.ex1$data1$clusters[2,]