[R] gamma parameter estimation [was "ks.test"]

[Martin Maechler]

>>>>> "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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._