### Compute conditional PDF using Kernel Density Estimators

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]
X=sesrain

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]

Y=nino3ses

#joint data
data=cbind(sesrain,nino3ses)

#where do you want to condition the PDF
# f(X | Y = Yc)  = f (X,Y=Yc) / f(Y = Yc)
# X is the seasonal rainfall
# Y is the seasonal NINO3
# Yc is the conditioning point..

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

#compute the joint PDF f(X,Y=Yc)
# for R version > 2.4 using the latest 'sm'
#get the joint PDF f(X,Y=Yc)
zz=sm.density(data,eval.points=xyc, display="none")

#get the marginal pdf at the conditioning point, i.e., Y = -0.5
# for R version > 2.4 using the latest 'sm'
fxc=sm.density(nino3ses,eval.points=Yc, display="none")

#get the conditinoal PDF  f(X,Y=Yc)/f(Y=Yc)
# for R version > 2.4 using the latest 'sm'
conpdf=diag(zz$estimate)/fxc$estimate
#conpdf=zz$estimate/fxc$estimate

plot(X1,conpdf,type="l",xlab="Rainfall ",ylab="Conditional PDF")

title(main=c("Conditional PDF - at ", Yc))