{

#To do a contigency table..
# Given aENSO forecast - compute the probabilities..

# read in the data file as a matrix of 13 columns
# first column is the year and the next 12 are the monthly
# rainfall values

rain=matrix(scan("aismr.txt"),ncol=13,byrow=T)

#pick the monthly rainfall values excluding the year.
rain1=rain[,2:13]
years=rain[,1]		#1871 - 1999		
N=length(years)		#of years
N1=N-1

#calculate monthly mean, variance and skews


zmean=colSums(rain1)/N

#or
zmean=apply(rain1,2,mean)
zvars=apply(rain1,2,var)
zskew=apply(rain1,2,skew)




#problem 5 - use boxplots for testing mean differences..
#read the NINO3 index..
test=matrix(scan("nino3-index.txt"),ncol=3,byrow=T)
nino3=matrix(test[,3],ncol=12,byrow=T)

#get the nino3 of the season
nino3ses=apply(nino3[,6:9],1,mean)

nino3years=unique(test[,1])		#1900 to 2004

#get the seasonal rainfall for the common years - i.e. 1900 - 1999
N=1999-1900+1

N2=1900-1871+1
N3=N2+N-1

sesrain=apply(rain1[,6:9],1,mean)
sesrain=sesrain[N2:N3]

datamat=cbind(sesrain,nino3ses[1:100])
xrain=c(min(sesrain)-1,quantile(sesrain,0.25),quantile(sesrain,0.75),max(sesrain)+1)
xnino=c(min(nino3ses)-1,-0.75,0.75,max(nino3ses)+1)
datamat=as.data.frame(datamat)
ztable=with(datamat, table(cut(datamat[,1],xrain, dig.lab=5), cut(datamat[,2],xnino,dig.lab=5)))

#the resulting table will have 3 x 3 = 9 boxes. The rows are Low, Mediaum
# and High rainfall categories and the columns are La Nina, Netural and
# El Nino categories.


#given ENSO forecast of 0.2,0.2,0.6	(La, Neut, El)
ensofcast=c(0.2,0.2,0.6)

#get the conditional probabilities
ztable1=ztable
ncat=3			#number of categories
for(i in 1:ncat){
ztable[,i]=ztable1[,i]/sum(ztable1[,i])
}

lflowc=ztable[1,]/sum(ztable[1,])
tpropl=sum(ensofcast*lflowc)

hflowc=ztable[3,]/sum(ztable[3,])
tproph=sum(ensofcast*lflowc)

#--------------------------------------------------------------------
#install library(gregmisc)
library(gplots)
library(akima)		#perhaps not required wth gplots..

i1=8
i2=9

 h2d <- hist2d(rain[,i1],rain[,i2],show=FALSE, same.scale=TRUE, nbins=c(20,30))

persp( h2d$x, h2d$y, h2d$counts,
   ticktype="detailed", theta=30, phi=30,
   expand=0.5, shade=0.5, col="cyan", ltheta=-30)

}