{
# 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

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

#plot the histogram
hist(x,xlab="Monthly rainfall", ylab="", probability=T, main="")
title(main=c(months[i],' rainfall'))
lines(xeval,xdensitynor)
lines(xeval,xdensitylnor,lwd=2)

}