[R] Maximum Likelihood Estimation

[Ben Bolker]

Todd Brauer <toddbrauer <at> yahoo.com> writes:

> 
> Using R, I would like to calculate algorithms to estimate coefficients á and â
within the gamma function:
> f(costij)=((costij)^á)*exp(â*costij).  I have its logarithmic diminishing line
data
> (Logarithmic Diminishing Line Data Table) and have installed R¢s Maximum
Likelihood Estimation
> package; however, I am unsure which method to apply in order to calculate the
algorisms (i.e.,
> Newton-Raphson Maximization, Nelder-Mead Maximization, etc.)  Any guidance you
all could provide
> would be appreciated.

  For a simple, well-behaved problem (like this one) it doesn't
really matter that much what algorithm you pick. optim()'s default
is Nelder-Mead (typically slightly more robust and less efficient).
Here's some sample code to do a least-squares fit:

x <- read.table("gammafit.txt")
names(x) <- c("x","y")
plot(x)
sum(x)
gfun <- function(p) {
    a <- p[1]
    b <- p[2]
    exp.y <- x$x^b*exp(-a*x$x)
    sum((exp.y-x$y)^2)
}
opt1 <- optim(gfun,par=c(a=0.1,b=1))
v <- opt1$par
curve(x^v["b"]*exp(-v["a"]*x),add=TRUE,col=2)

  However, if you really want to do MLE you may need to specify
your problem a little more carefully -- sum of squares will
find a reasonable answer but likelihood-based inference won't
be right.  Is this a survival problem (in which case you
might want to use the survival package)?  If you construct
your own likelihood function, you might need
to use pgamma(...,lower.tail=FALSE).

  Ben Bolker