to_add_to_pkg/PlotNLLbyM.r

In CCSBT/sbtr: sbtr

#######################################################
#Plots the likelihoods for different likelihood components and the total
#split out by different values of M10. 
#Idea and original code: Paige Eveson
#Code modified and rewritten by Trevor A Branch 
#  6 September 2009  tbranch@gmail.com
#Modified 24 June 2014 Trevor A Branch
#To split out by different values of M10 instead of steepness
#######################################################

#read in all the files into an object (takes time)
#source("get.all.data.r")
#data.M4constrained <- get.all.files(directory="M4constrained\\arc\\c1s1l1sqrt")

#get table of likelihood components and steepness values (quick)
#source("TableLikelihoodComponents.r")
#LikeM4const  <- likelihood.table(data.objects = data.M4constrained)
#LikeM4unconst <- likelihood.table(data.objects = data.M4unconst)
#LikeMpower <- likelihood.table(data.objects = data.Mpower)

#Plot the likelihoods by M10 values, using Paige's code and idea. 
plot.NLL.by.M <- function(NLLtable,caption) {
   par(mfrow=c(3,4),mai=c(0.1,0.2,0.2,0.2),omi=c(0.6,0.6,0.4,0))
   M10.vals <- sort(unique(NLLtable[,"M10"]))
   
   x <- as.vector(NLLtable[,"steepness"])
   x.lev <- sort(unique(x))
   ncolors <- length(x.lev)
   col.vec <- vector(length=length(x.lev))
   for (i in 1:length(x.lev)) {
      col.vec[x==x.lev[i]] <- i   
   }
   #print(col.vec)
   
   plot.nums <- c(1:9,11:13)
   nplots <- length(plot.nums)
   for(iii in 1:nplots) {  #the likelihood components
      if(iii==1) {  #do legend
         plot(x=1,y=1,type="n",axes=F,ann=F,ylim=c(0,2),xlim=c(0,2))
         legend(x=0,y=2.45,bty="n",horiz=T,
                pch=19,cex=1.3,col=(1:ncolors)+1,
                legend=paste("h =",x.lev),
                xpd=NA)
         par(new=T)
      }
      
      i <- plot.nums[iii]
      diffy <- max(NLLtable[,i]) - min(NLLtable[,i])
      
      if (diffy < 2.2) {
         midd <- mean(c(min(NLLtable[,i]),max(NLLtable[,i])))
         ylim <- c(midd-1.1,midd+1.1)
      }
      else {
         ylim <- c(min(NLLtable[,i]),max(NLLtable[,i])+0.1*diffy)
      }
      plot(x=jitter(NLLtable[,"M10"],amount=0.003),y=NLLtable[,i],xlab="",
                     ylab="",las=1,axes=F,ylim=ylim,pch=19,cex=0.8,col=col.vec+1)
      box()
      axis(2,las=1)
      if (iii>8) {
         axis(1,at=M10.vals)
      }
      coords <- par("usr")
      abline(h=mean(coords[3:4])-1,lty=2)
      abline(h=mean(coords[3:4])+1,lty=2)
      #par(xpd=T)
      #arrows(x0=coords[1],x1=coords[1],y0=coords[3],y1=min(coords[4],coords[3]+1),col="green",lwd=2,length=0)
      #par(xpd=F)
      mtext(outer=F,line=-2, colnames(NLLtable)[i],cex=1.2)
   }
   mtext(side=1,"M10",outer=T,line=3,cex=1.3)
   mtext(side=2,"Negative log likelihood",outer=T,line=2.5,cex=1.3)
   mtext(side=3,caption,outer=T,line=0.5,cex=1.3)
}