{

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

library(akima)

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

i1=8
i2=9
X=rain1[,i1]
Y=rain1[,i2]
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 x=max(X)/2

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

xc=max(X/2)
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



}