reports/genfig.R

   1 nvals <- c(5000,10000,100000,500000,1000000)
   2 for (link in c("logit","probit"))
   3 {
   4     for (d in c(2,5,10))
   5     {
   6         par(mfrow=c(2,2), lwd=2, cex.axis=2, cex=2)
   7         res <- list()
   8         for (n in nvals)
   9         {
  10             load(paste("multirun_",n,"_",d,"_",link,".RData",sep=""))
  11             res <- c(res, mr)
  12         }
  13         for (i in 1:2) #our, fm
  14         {
  15             for (j in 1:2) #beta, p,b
  16             {
  17                 if (j == 1) #beta
  18                 {
  19                     for (dim in 1:d)
  20                     {
  21                         ypts <- c()
  22                         for (k in 1:5)
  23                             ypts <- c(ypts, mean(sapply(1:length(res[[k]][[i]]), function(x) mean((res[[k]][[i]][[x]][2:(d+1),dim] - mr_params$beta[,dim])^2))))
  24                         plot(nvals, ypts)
  25                         if (dim < d)
  26                             par(new=TRUE)
  27                     }
  28                 }
  29                 else #p + b
  30                 {
  31                     for (rowidx in c(1,d+2))
  32                     {
  33                         ypts <- c()
  34                         ref <- if (rowidx==1) { mr_params$p } else { mr_params$b }
  35                         for (k in 1:5)
  36                             ypts <- c(ypts, mean(sapply(1:length(res[[k]][[i]]), function(x) mean((res[[k]][[i]][[x]][rowidx,] - ref)^2))))
  37                         plot(nvals, ypts)
  38                         if (rowidx==1)
  39                             par(new=TRUE)
  40                     }
  41                 }
  42             }
  43         }
  44     }
  45 }