#This code demonstrates the use of rjags to analyze loss data from Collins and Lowe 2005, See: http://sciencepolicy.colorado.edu/admin/publication_files/resource-2476-2008.02.pdf
#We examine total US losses due to hurricanes. 
#The predictors are Sept Sun Spot Number,  May June average North Atlantic Oscillation (NAO), August-October averaged Southern Oscillation Index (SOI), and August-October averaged North Atlantic Sea Surface Temperatures(SST)
#To use: 
#Put all code in the same directory
#set working directory to that directory.
#This code requires the files:
#jagsgpd.txt, jagsgpddata.txt, jagsgpdinit.txt,Covariates.txt,JagsExampleGpdfunctions.R ,GuiInterface.R, guilist.R
#This code demonstrates the use of jags.
#First we load in the data files into R, since jags uses R
#install.package("rjags" from Cran-R
#Function to read in the jags data and another to write out a list of data for use in command line jags. 
#The read is set to create environments and return the list, you can actually just return the environment and pass it to jags
#In order to use the .Rdata one requires that the code is set up for it.
#We have removed all code except the "jagsSampleMatrix" 
clean<-F
rerun<-F
if(exists("jagSampleMatrix")) reload<-F else reload<-T

if(rerun)
{
require(rjags)  #must have jags installed
source("JagsExampleGpdfunctions.R") 
jagsData<-read.jagsdata("jagsgpddata.txt")
jagsInits<-read.jagsdata("jagsgpdinit.txt")
testCovariates<-data.frame(read.jagsdata("Covariates.txt"))
#Create a list of the random number generators. We will use all of them in our test
RNGs<-c("base::Wichmann-Hill","base::Marsaglia-Multicarry","base::Super-Duper","base::Mersenne-Twister") 
options(jags.pb="gui") #sets up a gui for the progress bar
#Initialize Jags Model
jagsmodel<-jags.model(file="jagsgpd.txt",data=jagsData,inits=lapply(RNGs,function(x) {xx<-jagsInits; xx$.RNG.name=x ; xx}),n.chains=length(RNGs),n.adapt=4000 )#Sets up the model call
#Generate samples (could use 25,000 but we can get 1000 easily hear: We sample 4 chains with same initializations but different random number generators.
jagSamples<-coda.samples(jagsmodel,c("Rate","Sigma","Xi","deviance"),n.iter=25000) #change to 25000
##Since we used arrays we wanted to change the names of the variables so here are the actual predictors in the model
parmaskgpd=list(LogRate=c("INT","sst","nao","soi","ssn"),LogSigma=c("INT","sst","nao"),Xi=c("INT","soi")); parnamesgpd<-c(unlist.names(parmaskgpd),"deviance")#Names chosen to match sorted results from coda.samples...
#Here we attach new names to the model regression coeffficients which were Rate[1], Rate[2], Sigma[1], Sigma[2],Sigma[3], Xi[1], Xi[2], deviance to
#LogRate.INT , LogRate.soi , Sigma.INT , Sigma.nao , Sigma.sst , Xi.INT , Xi.soi , deviance 
jagSamples<-lapply(jagSamples,function(x){ colnames(x)<-parnamesgpd; x}) ; class(jagSamples)<-"mcmc.list"
print(jagStats<-summary(jagSamples))
#Plot the densities, I created a version of plot.mcmc to manage my own output size and to perform a density plot on all the data
jplots(tplot,jagSamples,trace=F,title="Posterior Parameter densities for Collins and Lowe loss data using GPD model",mfrow=c(3,4))
#plot the BGR diagnostics
jplots(gelman.plot,jagSamples,title="BGR convergence statistics for Collins and Lowe loss data using GPD model",mfrow=c(3,4)) 
#plot the autocorrelation with maximum lag=100 for the first sample
jplots(autocorr.plot,jagSamples[[1]],paste("Parameter sample Autocorrelation plots forCollins and Lowe loss data using GPD model"),lag.max=100,mfrow=c(3,4))
jagSampleMatrix<-as.matrix(jagSamples) #also as.array to separate each set of samples Merges chains for each variable
#The following is an example of how you might use the matrix in predictions:
testCovariates<-data.frame(read.jagsdata("Covariates.txt"))

#CL2005 is Collins and Lowe 2005 
if(clean)
{ 
saveobjects<-c("jagSampleMatrix","jagsData","parmaskgpd","testCovariates","tomboxes","plotextreme","jagsgpdposterior","samplergpd","maketext","returnlevels.sample","as.array.list","unlist.names","make.add.text")
save(list=saveobjects,file="SavedObjects.R")
rm(list=ls())
reload<-T
}

}
if(reload) load("SavedObjects.R")
jagsPosteriorSamples<-jagsgpdposterior(x=jagSampleMatrix,thin=10, newcov=testCovariates[c(2,3,4,6),],N=NULL,u=jagsData$u,mask=parmaskgpd,yp=c(1,2,5,10,20,50,100,200,500,1000),mfrow=c(2,2),yrange=c(7,16),yl=18,type="CL2005",etafun=c(exp,exp,identity),plot="Dot")
#
#use testCovariats to generate mu and sd, inputs are based on mean and standard deviations:
musd<-rbind(mu=testCovariates[3,],sd=abs(testCovariates[2,]-testCovariates[3,]))
require(gWidgets)
source("GuiInterface.R")
source("guilist.R")
guiToolkit("RGtk2") #It limps in tcl/tk but it will work. guiToolkit("tcltk")
require(cairoDevice)
options(guiToolkit="RGtk2")
guiwindow<-gwindow("Posterior distribution of common logrithm of loss in 2005 Dollars",width=300,height=500)
guiCommandLine<-gcommandline(console=T)
tomgeneric(lst=guilist,cli=guiCommandLine,container=guiwindow,toolkit=guiToolkit())
posteriorgraph<-ggraphics(container=guiwindow,ps=12,width=600,height=400)