This preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Unformatted text preview: # calculation of MLE for Gamma distribution # illinois 1960 rainfall data (see chapter 10) # This function computes the method of moments estimator and its approximate variance 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.2mu.1^2 V.hat = var(cbind(x,x^2)) lam.hat = mu.1/(mu.2mu.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.2mu.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.2mu.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) } # the MLE needs to be computed using an interative method, as there is no explicit solution to the MLE # it appears easier to use the minimizer optim, but two R functions to find the MLE are given here # this function gives the negative log likelihood in a form to be used in the nonlinear minimier nlm NLik = function(a,x){ lam = a[1] alpha = a[2] n = length(x) L = n*alpha*log(lam) + (alpha1)*sum(log(x))  lam*sum(x)  n*log(gamma(alpha)) L = (1)*L # .grad is a 1 by 2 matrix that is the derivatives of L .grad = array(0, c(1,2), list(NULL, c("lam", "alpha")))....
View
Full Document
 Spring '11
 qqqq

Click to edit the document details