#fit a Bivariate normal PDF and conditional..

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

library(akima)
library(mvtnorm)

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("http://civil.colorado.edu/~balajir/r-session-files/nino3-index.txt"), ncol=3, byrow=T)

#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=sesrain
Y=nino3ses
data=cbind(X,Y)

N = length(X)

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

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

rho=cor(X,Y)

nx=60
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)

k=0
for(j in 1:nx){

	for(i in 1:ny){
k=k+1
xyeval[k,1]=xeval[j]
xyeval[k,2]=yeval[i]
nfact=1/(2*pi*sigx*sigy*sqrt(1-rho^2))
xs=(xyeval[k,1]-xm)/sigx
ys=(xyeval[k,2]-ym)/sigy

jpdf[k]=exp(-(xs*xs + ys*ys - 2*rho*xs*ys)/(2*(1-rho^2)))
jpdf[k]=jpdf[k]*nfact

}

}

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

persp(zz,ticktype="detailed", theta=30, phi=30, expand=0.5, shade=0.5, col="cyan", ltheta=-30)

#or
contour(zz,xlab="X", ylab="Y")
points(X,Y)


xyeval = expand.grid(x=xeval, y=yeval)
jpdf = dmvnorm(xyeval, mean=colMeans(data), sigma=var(data))
zz=interp(xyeval[,1],xyeval[,2],jpdf)
persp(zz,ticktype="detailed", theta=30, phi=30, expand=0.5, shade=0.5, col="cyan", ltheta=-30)

#or
contour(zz,xlab="X",yalb="Y")
points(X,Y)


#conditional slize at Y = Yc

#f(x | y) = f(x,y=Yc)/f(y=Yc)

Yc=-0.5
Yc=0.5
X1 = sort(sesrain)
#xyc=cbind(X1, rep(Yc,n))

#or to get a finer grid..
xyc=cbind(xeval, rep(Yc,nx))


fxyc = dmvnorm(xyc, mean=colMeans(data), sigma=var(data))

fy = dnorm(Yc,mean=mean(Y), sd=sd(Y))		#PDF of Y evaluated at Yc

fxcond = fxyc/fy
plot(xeval, fxcond, type="l", xlab="Rainfall", ylab="conditional PDF")

#nfact=1/(sqrt(2*pi)*sigx)
#fxc=nfact*exp(-(xc-xm)*(xc-xm)/(2*sigx*sigx))
#fyc=1:ny
#for(i in 1:ny){
#nfact=1/(sqrt(2*pi)*sigy)
#fyc[i]=nfact*exp(-(yeval[i]-ym)*(yeval[i]-ym)/(2*sigy*sigy))}
#fyc=fyc/fxc
#
#plot(yeval,fyc,type="l",xlab="Rainfall",ylab="conditional PDF")