###################################################
### chunk number 1: startchunk
###################################################
require(affy)
require(affydata)
data(Dilution)


###################################################
### chunk number 2: RPA.pointestimate
###################################################
require(RPA)
sets <- geneNames(Dilution)[1:5]
rpa.results <- RPA.pointestimate(Dilution,sets,cind=1)


###################################################
### chunk number 3: rpa.res eval=FALSE
###################################################
## rpa.results <- RPA.pointestimate(Dilution,cind=1)


###################################################
### chunk number 4: visu
###################################################
set <- rownames(rpa.results$d)[[1]]
d <- rpa.results[[set]]$d
s2 <- rpa.results[[set]]$sigma2
plot(rpa.results[set,])


###################################################
### chunk number 5: barplots
###################################################
plot(rpa.results[set,])


###################################################
### chunk number 6: Smat
###################################################
Smat <- RPA.preprocess(Dilution,cind=1)


###################################################
### chunk number 7: 
###################################################
set = "1000_at"
S <- t(Smat$fcmat[pmindex(Dilution, set)[[1]], ])


###################################################
### chunk number 8: 
###################################################
res <- RPA.iteration(S, sigma2.guess = rep(1,ncol(S)))


###################################################
### chunk number 9: 
###################################################
my.priors = initialize.priors(Dilution, sets, alpha = 1e-6, beta = 1e-6)


###################################################
### chunk number 10: 
###################################################
set = "1000_at"
probe.idx = 5 # set priors for this probe within the probeset
my.priors[[set]]$beta[[probe.idx]]  = 10
my.priors[[set]]$alpha[[probe.idx]] = 10


###################################################
### chunk number 11: 
###################################################
rpa.results <- RPA.pointestimate(Dilution,sets,
                                 priors=my.priors,
                                 sigma2.method = "basic",
                                 d.method = "basic")


###################################################
### chunk number 12: 
###################################################
barplot(rpa.results[[set]]$sigma2)


###################################################
### chunk number 13: 
###################################################
eset = rpa2eset(rpa.results)


###################################################
### chunk number 14: 
###################################################
set.seed(24)
d.true <- rnorm(300,mean=0,sd=4)
sigma2.true <- c(1,1,1,1,2,2,3,3,3,4,6)
S <- NULL
for (s2 in sigma2.true) {S <- cbind(S,rnorm(length(d.true),mean=d.true,sd=sqrt(s2)))}


###################################################
### chunk number 15: 
###################################################
set.seed(234)
res <- RPA.iteration(S, sigma2.guess = rep(1,ncol(S)))


###################################################
### chunk number 16:  eval=FALSE
###################################################
## par(mfrow=c(2,2))
## plot(res$d,d.true,main="True vs. estimated d")
## plot(res$sigma2,sigma2.true, ,main="True vs. estimated sigma2")