[R] gamma parameter estimation [was "ks.test"]
Martin Maechler
maechler at stat.math.ethz.ch
Mon Nov 12 18:49:34 CET 2001
>>>>> "AOlinto" == aolinto <aolinto at bignet.com.br> writes:
AOlinto> Dear Dr. Maechler, Thanks for your e-mail.
MM> As the function name suggests and help(gamma.shape) clearly says,
MM> its estimating the shape parameter alpha, and *ITS* standard error
MM> which is not at all the same as the scale parameter you need (and
MM> in your case get from coef(Lt.fit) !
{to R-help; I allow myself to CC this back to R-help again..}
AOlinto> As you stated, I was not using the scale parameter, but the SE
AOlinto> of the shape parameter.
AOlinto> Nevertheless, in my exemple, using coef(Lt.fit) I got the
AOlinto> value 0.01703055(Intercept) which doesn't seem to be the scale
AOlinto> parameter either. I read carefully some books and the R help
AOlinto> but I couldn't find a solution (probably because I'm not a
AOlinto> statistician and, as biologist, I often have problems with
AOlinto> mathematical language - but I'm trying to do my best).
AOlinto> Would you please explain me how to get the scale parameter
AOlinto> from coef?
Here is a valid R script I did for you [ ~/R/MM/STATISTICS/gamma-mle.R ]:
the comments should contain all the explanations:
Note that it also contains a very simple ("method of moments") procedure
for estimation of gamma parameters. `Method of moments' is sub-optimal
theoretically, but in this case has the advantage to be much easier.
----------------------
### Create reproducible sample data (large N such that estim. ~= true):
### try a smaller N, e.g., N = 100, to see that MLE is better
set.seed(35)
hist(xx <- rgamma(10000, shape = 2, scale = 0.1))
## =alpha =sigma = 1/rate
## help(rgamma): shape = alpha =: a, scale = sigma =: s
## -------------> E[X] = a * s; Var[X] = a * s^2
c(mean = (m <- mean(xx)), var = (v <- var(xx)))
## 1) Method of moment estimator :
## =============================
sigma0 <- v/m
alpha0 <- m / sigma0
c(alpha0 = alpha0, sigma0 = sigma0)
##- alpha0 sigma0
##- 1.9794580 0.1001518
all(alpha0*sigma0 ^ (1:2) == c(m,v)) # TRUE : moments match
## 2) Better method: MLE -- via GLM and gamma.shape():
## ------------------------------------------------
## What does the Gamma family do?
(aicG <- Gamma()$aic) # is -2 * logLik(.) + 2
## i.e.,
## logLik() = sum(dgamma(y, 1/disp, scale = mu * disp, log = TRUE) * wt)
## or the *unweighted* likelihood
## Lik() = prod( dgamma(y, shape = 1/disp, scale = mu * disp) )
Gamma()$link ## "inverse", i.e. the intercept below is c0 = 1/mu
## c0 = 1 / mu = 1/ (sigma * alpha) = "rate" / alpha
## or equivalently : sigma = 1/(c0 * alpha)
gx <- glm(xx ~ 1, family = Gamma)
(sgx <- summary(gx))# see also "corrected" summary below
(cc <- coef(gx))# 5.044 +/- 0.036
ccV <- sgx $ coefficients["(Intercept)" , ]
cc == ccV[1] # TRUE
ccV["Estimate"] + c(-2,2)* ccV["Std. Error"] # [ 4.97, 5.12 ]
library(MASS)
(gs <- gamma.shape(gx, eps = 1e-7, verbose=TRUE)) # alpha = 1.974 (+/- 0.026)
gs$alpha + c(-2,2)* gs$SE # [1.922, 2.025 ]
## corrected summary using mle dispersion (but here, change is minuscule!):
(sgxC <- summary(gx, dispersion = 1 / gs$alpha))
## = summary(gx, dispersion = gamma.dispersion(gx))
ccV. <- sgxC $ coefficients["(Intercept)" , ]
ccV - ccV. ## no difference in this case
## `MLE' sigma:
(sigma <- 1 / (coef(gx) * gs$alpha))# = 0.1004 in the example
AOlinto> Thanks for your attention,
AOlinto> Antonio Olinto
you're welcome!
Martin Maechler <maechler at stat.math.ethz.ch> http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum LEO D10 Leonhardstr. 27
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
phone: x-41-1-632-3408 fax: ...-1228 <><
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list