We will follow the series of examples in chapter 8 which use the data set ILLRAIN from 
the text's diskette.
In a homework problem (chapter 8, problem 46) you are asked to perform the same 
calculations using the data set
whales.txt.
The dataset ILLRAIN.txt is available on my webpage (as is whales.txt). 

Rainfall of summer storms was measured by a network of rain gauges 
in southern Illinois for the years 1960-1964.

The dataset as listed in the file uses the symbol "*" to denote missing data.
The "R" package does not like that so we take care of it with the option na.strings="*".

Notice also, that I allowed "R" to choose my header names for me.

illrain<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/ILLRAIN.txt",header=F,na.strings="*")

or if the dataset is in your working directory:


illrain<-read.table("U:/Math 152/ILLRAIN.txt",header=F,na.strings="*")

(Here is  the corresponding line that I used on my machine at home)

illrain<-read.table("C:/Documents and Settings/Mike/Desktop/Math 152/ILLRAIN.txt",header=F,na.strings="*")

attach(illrain)

rain<-c(V1,V2,V3,V4,V5)


As it stands, the vector "rain" contains many entries of the form "NA".
Though it may not be necessary for what follows, they make me uncomfortable, so I will remove them.
The logical negation operator is "!" and "is.na(rain)" is a list of the values TRUE and FALSE with TRUE in 
the positions occupied by "NA".

rain<-rain[!is.na(rain)]

rain

hist(rain)

The histogram is skewed, indicating that a Gamma distribution might be a good choice 
for a fit to the data. (see pg. 250)

We will first fit a Gamma distribution to the data by the method of moments.


varhat<-mean(rain^2)-mean(rain)^2

lambdahat<-mean(rain)/varhat

alphahat<-lambdahat* mean(rain)

f<-function(x){((lambdahat^alphahat)/gamma(alphahat))*x^(alphahat-1)*exp(-lambdahat*x)}

hist(rain,probability=T)

curve(f(x),add=T)


(alternately one can use something like

curve(dgamma(x,shape=5,rate=1),xlim=c(1,20),ylim=c(0,.2))

)


The fit is reasonably good in this case (but need not be in general). 
------------------------------------------------------------------------------------
Bootstrapping with  the method of moments:

Recall for the method of moments

varhat<-mean(rain^2)-mean(rain)^2

lambdahat<-mean(rain)/varhat

alphahat<-lambdahat* mean(rain)

alphahat

lambdahat


resamp<-function(x){rgd<-rgamma(227,shape=.3779,rate=1.6842)
                    mean(rgd^2)-mean(rgd)^2->var
                    lambdastar<-mean(rgd)/var
                    alphastar<-lambdastar*mean(rgd)
                    lambdastar
 }

hist(sapply(1:1000,resamp))


or,


replicate(1000,

{rgd<-rgamma(227,shape=.44079,rate=1.964374)
                    mean(rgd^2)-mean(rgd)^2->var
                    lambdastar<-mean(rgd)/var
                    alphastar<-lambdastar*mean(rgd)
                    lambdastar
 }
)


I personally find the last form to be more cumbersome so I will 
stick with the functional form.

data<-sapply(1:1000,resamp)
slambdahat<-sqrt((1/1000)*sum((data-mean(data))^2))
slambdahat

Now for the alphahat calculation:

resamp<-function(x){rgd<-rgamma(227,shape=.3779,rate=1.6842)
                    mean(rgd^2)-mean(rgd)^2->var
                    lambdastar<-mean(rgd)/var
                    alphastar<-lambdastar*mean(rgd)
                    alphastar
 }

data<-sapply(1:1000,resamp)
salphahat<-sqrt((1/1000)*sum((data-mean(data))^2))
salphahat


These last values should be compared with those in the text and with different 
determinations.
-----------------------------------------------------------------------------------------