#fit a Bivariate normal PDF and conditional..

#install library (gregmisc)
#then
#library(gplots)

library(akima)
library(sm)

rain=matrix(scan("http://civil.colorado.edu/~balajir/r-session-files/aismr.txt"),ncol=13,byrow=T)


#rain=matrix(scan("aismr.txt"),ncol=13,byrow=T)
N=1999-1900+1

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

sesrain=apply(rain[,7:10],1,mean)
sesrain=sesrain[N2:N3]

n=length(sesrain)

#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 summer (June - Sep) season
nino3ses=apply(nino3[,6:9],1,mean)

nino3ses=nino3ses[1:n]

X=nino3ses
Y=sesrain

X=sesrain
Y=nino3ses

xm=mean(X)
sigx=sd(X)

ym=mean(Y)
sigy=sd(Y)

rho=cor(X,Y)

nx=50
ny=50

nxy=nx*ny
xeval=seq(min(X),max(X),length=nx)
yeval=seq(min(Y),max(Y),length=ny)

jpdf=1:nxy
xyeval=matrix(0,nrow=nxy,ncol=2)
xyeval=expand.grid(x=xeval, y=yeval)

data=cbind(X,Y)
jpdf = sm.density(data, eval.points=xyeval, eval.grid=FALSE,
display="none")

zz=interp(xyeval[,1],xyeval[,2],jpdf$estimate)

persp(zz,ticktype="detailed", theta=30, phi=30, expand=0.5, shade=0.5,
col="cyan", ltheta=-30,xlab="X", ylab="Y", zlab="PDF")