library(MASS)
library(fitdistr)
library(sm)


# read in the data file as a matrix of 13 columns
# first column is the year and the next 12 are the monthly
# rainfall values

#Fit, Normal, Lognormal, Gamma, Weibull and also Nonparametric density
# estimates to the data..

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

#calculate monthly mean, variance and skews

zmean=apply(rain1,2,mean)
zvars=apply(rain1,2,var)
zsd=apply(rain1,2,sd)
zskew=apply(rain1,2,skew)

#fit a Normal distribution for each monthly rainfall and plot them on top of the 
#histogram

months=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")

par(mfrow=c(3,4))

for(i in 1:12){
x=rain1[,i]
xeval=seq(max(0,min(x)-zsd[i]),max(x)+zsd[i],length=100)
fact=1/(sqrt(2*pi)*zsd[i])
xdensitynor=fact*exp(-(xeval-zmean[i])*(xeval-zmean[i])/(2*zvars[i]))

#or

xdensitynor=dnorm(xeval,mean=zmean[i],sd=zsd[i])
xdensitylnor=dlnorm(xeval,meanlog=mean(log(x)), sdlog=sd(log(x)))

#Gamma
zz=fitdistr(x,"gamma")

xdensitygamma=dgamma(xeval,shape=zz$estimate[1],scale=1/zz$estimate[2])

zz=fitdistr(x,"weibull")
xdensityweibull=dweibull(xeval,shape=zz$estimate[1],scale=zz$estimate[2])
 
#plot the histogram
zz=hist(x,freq=FALSE,plot=FALSE)
hist(x,xlab="Monthly rainfall", ylab="", freq=FALSE,
main="",ylim=range(c(zz$density,xdensitynor, xdensitylnor, xdensitygamma, xdensityweibull)))
title(main=c(months[i],' rainfall'))
lines(xeval,xdensitynor,col="red")
lines(xeval,xdensitylnor,col="orange")
lines(xeval,xdensitygamma,col="blue")
lines(xeval,xdensityweibull,col="green")

#Nonparametric Kernel density estimation..

lines(xeval, sm.density(x,eval.points=xeval,add=TRUE)$estimate, lwd=4)
box()

}