ALLPLOT <- T
### T for all 2^5*3^2 plots
### F for six randomly selected plots

source("analyze.q")

### HYPERPARAMETERS
sd.alpha <- c(10,10)
sd.beta <-  c(10,10,0.5)
shp.lam <- 0.01; scl.lam <- 0.01
shp.sig <- 0.01; scl.sig <- 0.01
shp.tau <- 0.01; scl.tau <- 0.01

shp.tau <- 5; scl.tau <- 1.25

### DESIGN SETTINGS

VARRAT.VAL   <- c(0.25, 0.5)
RHO.XS.VAL   <- c(0.25, 0.5, 0.75)
RHO.YX.S.VAL <- c(0.15, 0.5)
RHO.YS.X.VAL <- c(0, 0.15, 0.3)

if (!ALLPLOT) {
  postscript("/usr/local/data/gustaf/instrument/scriptA.ps")
  par(mfrow=c(3,4)) }

### j1,j2,j3,j4 index the two-level factors

set.seed(133); set.seed(17)
RANDPLOT <- sample( c(rep(F, 2^4*3^2 - 8), rep(T, 8) ) )
tmp <- RANDPLOT; tmp[tmp==T] <- c(rep(F,4),rep(T,4))
RANDPLOT <- c(RANDPLOT,tmp)
plotndx <- 0

for (j1 in 1:2) {  for (j2 in 1:2) { for (j3 in 1:2) { for (j4 in 1:2) {
for (j5 in 1:2) {
  
  ### which prior on beta3
  if (j1==1) { sd.beta[3] <- 0.5 }
  if (j1==2) { sd.beta[3] <- 5   }
  
  ### which prior on tau2
  if (j2==1) { shp.tau <- 4; scl.tau <- 1}
  if (j2==2) { shp.tau <- .01; scl.tau <- 0.01 }

  ### how much meas err
  tau2 <- VARRAT.VAL[j3]

  ### sign of beta[3]
  sgn3 <- c(-1,1)[j4]

  ### size of target effect
  RHO.YX.S <- RHO.YX.S.VAL[j5]

  if (ALLPLOT) {
    postscript(paste("/usr/local/data/gustaf/instrument/",
                     c("L","H")[j1],
                     c("L","H")[j2],
                     c("L","H")[j3],
                     c("L","H")[j5],
                     c("L","H")[j4],
                     ".ps",sep=""))
    par(mfrow=c(3,3)) }

  for (i1 in 1:3) { for (i2 in 1:3) {

  RHO.XS <- RHO.XS.VAL[i1]
  RHO.YS.X <- RHO.YS.X.VAL[i2]

  ### three innocous choices
  beta  <- c(0, NA, NA)
  alpha <- c(0, NA)
  sig2 <- 1

  ### remaining design choices

  alpha[2] <- RHO.XS; lam2 <- 1-RHO.XS^2

  beta[3] <- sgn3*sqrt ( sig2 / (1-RHO.XS^2) * RHO.YS.X / (1-RHO.YS.X) )

  beta[2] <- sqrt ( sig2 / (1-RHO.XS^2) * RHO.YX.S / (1-RHO.YX.S) )

  tmp <- analyze(alpha, beta, tau2, lam2, sig2,
                 sd.alpha, sd.beta,
                 scl.tau, shp.tau, scl.lam, shp.lam, scl.sig, shp.sig,
                 GRLEN=2500)

  plotndx <- plotndx+1
  if (ALLPLOT || RANDPLOT[plotndx]) {

    if (ALLPLOT)  { yl <- max(tmp$pstdens) }
    if (!ALLPLOT) { yl <- 25               }
    plot(-100,-100, xlim=c(-1.5,1.5), 
               ylim=c(-0.22,1.05)*yl, col="blue",
               xlab=expression(beta[3]), ylab="Density")

  title(paste(c("L","H")[j1],
             c("L","H")[j2],
             c("L","H")[j3],
             c("L","M","H")[i1],
             c("L","H")[j5],
             c("L","H")[j4],
             c("L","M","H")[i2], sep=""))

 # if (tmp$lft>((-3)*sd.beta[3])) {
 #   rect(-3*sd.beta[3],-0.22*yl,tmp$lft,-0.12*yl,
 #        col=gray(.5),border=F) }

 # if (tmp$rht<((3)*sd.beta[3])) {
 #   rect(tmp$rht,-0.22*yl,3*sd.beta[3],-0.12*yl,
 #        col=gray(.5),border=F) }

  rect(tmp$lft, -0.22*yl, tmp$rht, -0.12*yl, col=gray(0.5), border=F)
    
  gr2 <- seq(from=(-3)*sd.beta[3], to=3*sd.beta[3], length=500)
  points(gr2, dnorm(gr2, sd=sd.beta[3]),type="l", col=gray(.5))

  points(tmp$gr, tmp$pstdens, type="l")

  points(beta[3], -0.08*yl, pch=17, cex=1.35)   }


} }

if (ALLPLOT) {  
graphics.off()  }

} } } } }

if (!ALLPLOT) {
  graphics.off() }