# sandm77.r                  November 2007, 2008
 #
 # replicating Swallow & Monahan simulation study
 #
 # read in function code for 'estims'
  source("estims.r")
 #
 # replication pattern
  ni <- c(3,3,5,5,7,7)      # this is P4
  a <- length(ni)
  N <- sum(ni)
 #
 # set seed
   set.seed(5151917 + 4 )   # for P4
 #
 # 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)
   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))
   c(nrow(ests),ncol(ests))
 # save results in file for further analysis
   write(ests,"estsP47.dat",ncolumns=9)       # goofy tranposition!
   warnings()
 #
 # print a few to see what it looks like
   print(t(ests)[1:5,])
   print(t(ests)[695:700,])
 # some quick diagnostics
   table(ests[7,],ests[8,])   # ml vs reml convergence
   table(ests[7,],ests[9,])   # ml convergence vs ratio
   table(ests[8,],ests[9,])   # reml convergence vs ratio
 # done
  rm(list=ls())
  q()