# sandm10.r                  November 2007, 2008
 #
 # replicating Swallow & Monahan simulation study
 #  -- this time, the whole thing
 #
 # read in function code for 'estims'
  source("estims.r")
 #
 # replication patterns
  nilist <- list( c(3,5,7), c(1,5,9), c(3,3,5,5,7,7), c(1,1,5,5,9,9),
         c(1,1,1,1,13,13), c(3,3,3,5,5,5,7,7,7), c(1,1,1,5,5,5,9,9,9),
         c(1,1,1,1,1,1,1,19,19), c(2,10,18), c(3,15,27) )
 #
 # initialize estall
   estall <- matrix(0,10,0)
 # set seed
   set.seed(5151917)          # usual seed
 #
 #                            main loop start
  for (pattern in (1:10)) {
 # 
  ni <- nilist[[pattern]]
  a <- length(ni)
  N <- sum(ni)
  print(ni)
 #
 # generate data
   nrep <- 100      # only a few replicates
   sigmaa <- rep(c(0,1,2,5,10,20,50)/10,rep(nrep,7)) # set true values
   sigmae <- 1      # of course
 #
 # first mean vectors
   ybars <- matrix(rnorm(a*nrep*7),a,nrep*7)      # 7 levels of sigmaa
   ybars <- ybars*sqrt(rep(sigmaa,a)+sigmae/ni)
 #
 # now sse/1 has chi-square distribution
   sses <- rchisq(nrep*7,(N-a)) # with (N-a) df
 #
 # put ybar and sse together as single input vector
 #   but columns of matrix  
   data <- rbind(ybars,unname(sses))    # drop labels
 # 
 # apply function to cols of matrix
   ests <- apply(data,2,estims)
   ests <- rbind(ests,unname(sigmaa),rep(pattern,ncol(ests)))
   print(pattern,nrow(ests),ncol(ests))
   estall <- cbind(estall,ests)
                          }   # end of pattern loop       
 #
 # save results in file for further analysis
   write(estall,"estsall.dat",ncolumns=10)       # goofy tranposition!
 #
 # print a few to see what it looks like
   print(t(estall)[1:5,])
 # some quick diagnostics
   table(estall[7,],estall[8,])   # ml vs reml convergence
   table(estall[7,],estall[9,])   # ml convergence vs ratio
   table(estall[8,],estall[9,])   # reml convergence vs ratio
   table(estall[7,],estall[10,])  # ml convergence vs pattern
   table(estall[8,],estall[10,])  # reml convergence vs pattern
 # done
  rm(list=ls())
  q()