{

#fit a Bivariate normal PDF and conditional..

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

library(akima)


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
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)

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

	for(j in 1:ny){
k=k+1
xyeval[k,1]=xeval[i]
xyeval[k,2]=yeval[j]
nfact=1/(2*pi*sigx*sigx*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)))


}

}

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)


#conditional slize at xc

xc=-0.5

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

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")


}