2.2 [R] Fit a Gamma distribution to Illinois rainfall data (Jan7)

2.2 [R] Fit a Gamma distribution to Illinois rainfall data (Jan7)

Info iconThis preview shows pages 1–2. Sign up to view the full content.

View Full Document Right Arrow Icon
# January 4, 2010 # calculation of MLE for Gamma distribution # illinois 1960 rainfall data (see chapter 10) # Method of moments fit ; discuss later gamma.moment.est.2 = function (x) { note = "Fit a Gamma distribution by method of moments" mu.1 = mean(x) mu.2 = mean(x^2) var.x = mu.2-mu.1^2 V.hat = var(cbind(x,x^2)) lam.hat = mu.1/(mu.2-mu.1^2) alpha.hat = mu.1*lam.hat # estimate of variance for lam.hat based on delta method a2 = array( c(mu.1^2+mu.2, -mu.1)/((mu.2-mu.1^2)^2) ,c(2,1)) tau.lam = t(a2)%*%V.hat%*%a2 a1 = array( c( 2*mu.1*mu.2, -mu.1^2)/((mu.2-mu.1^2)^2) ,c(2,1)) tau.alpha = t(a1)%*%V.hat%*%a1 list(theta.hat=c(lam.hat,alpha.hat),mu.1=mu.1,mu.2=mu.2,V.hat=V.hat,tau1=tau.alpha,tau2=tau.lam) } # read in the data rain.dat = scan("illinois60.txt") # fit parameters using method of moments rain.mom.fit = gamma.moment.est.2(rain.dat) # point estimate using method of moments theta.tilde = rain.mom.fit$theta.hat hist(rain.dat , freq = F) # not a very useful plot; need different break points
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: breaks = seq(0, 2.2 , .1) par(oma = c(0,0,4,0)) hist(rain.dat , freq = F, breaks = breaks, xlab = "rainfall (inches)", main = "Illinois Rainfall Histogram") # make a gamma density overlay on the histogram alpha.tilde = theta.tilde[2] lam.tilde = theta.tilde[1] x.pts = seq(0.01, 2.2, .01) d.gam.fit = dgamma( x.pts , shape = alpha.tilde , rate = 1/lam.tilde) lines( x.pts , d.gam.fit, lty = 2) # add some text in the margin on side 3 mtext("Gamma fit overlay", side = 3) # rainfall quantiles based on fitted model pr = c(.5 , .75 , .9 , .95 , .99) q.rain = qgamma( pr , shape = alpha.tilde , rate = 1/lam.tilde) rain.pred = cbind( pr , q.rain) print( rain.pred , digits = 3) Page 1 of 2 25/01/2011 http://www.stats.uwo.ca/faculty/kulperger/Stat3858/Computing/RScripts/Illinois-rain-1.txt Page 2 of 2 25/01/2011 http://www.stats.uwo.ca/faculty/kulperger/Stat3858/Computing/RScripts/Illinois-rain-1.txt...
View Full Document

Page1 / 2

2.2 [R] Fit a Gamma distribution to Illinois rainfall data (Jan7)

This preview shows document pages 1 - 2. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online