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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
