You asked for a hint.
> library(MASS)
> x <- c(20.0, 23.9, 20.9, 23.8, 25.0, 24.0,
21.7, 23.8, 22.8, 23.1, 23.1, 23.5,
23.0, 23.0)
> fitdistr(x, "gamma")
shape rate
316.563872 13.780766
(119.585534) ( 5.209952)
To do it with more general and elementary tools, look at ?optim
nls(...)? Not relevant at all, no matter how it feels.
Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary): +61 7 3826 7304
Mobile: +61 4 8819 4402
Home Phone: +61 7 3286 7700
mailto:[EMAIL PROTECTED]
http://www.cmis.csiro.au/bill.venables/
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
On Behalf Of Mohammad Ehsanul Karim
Sent: Sunday, 27 January 2008 4:48 PM
To: r-help@r-project.org
Subject: [R] Likelihood optimization numerically
Dear List,
I am not sure how should i optimize a log-likelihood numerically:
Here is a Text book example from Statistical Inference by George
Casella, 2nd
Edition Casella and Berger, Roger L. Berger (2002, pp. 355, ex. 7.4 #
7.2.b):
data = x = c(20.0, 23.9, 20.9, 23.8, 25.0, 24.0, 21.7, 23.8, 22.8,
23.1, 23.1, 23.5, 23.0, 23.0)
n <- length(x)
# likelihood from a 2 parameter Gamma(alpha, beta), both unknown
llk = -n*log(gamma(alpha)) - n*alpha*log(beta) + (alpha -
1)*(sum(log(x))) - (sum(x))/beta
# analytic 1st derivative solution w.r.t alpha, assuming beta known
# by putting MLE of beta = sum(x)/(n*alpha)
# (to simplify as far as possible analytically)
llk.1st = - n*digamma(alpha) -n*(log(sum(x)/(n*alpha))+1) +
(sum(log(x)))
It feels like i should use
nls(... , trace=T, start=c(alpha=...),nls.control(maxiter=100,tol=.1))
but not sure "how".
Can anyone provide me hint?
Thank you for your time.
Ehsan
http://www.youtube.com/profile_play_list?user=wildsc0p
> R.Version()
$platform
[1] "i386-pc-mingw32"
$arch
[1] "i386"
$os
[1] "mingw32"
$system
[1] "i386, mingw32"
$status
[1] ""
$major
[1] "2"
$minor
[1] "6.1"
$year
[1] "2007"
$month
[1] "11"
$day
[1] "26"
$`svn rev`
[1] "43537"
$language
[1] "R"
$version.string
[1] "R version 2.6.1 (2007-11-26)"
________________________________________________________________________
____________
Be a better friend, newshound, and
______________________________________________
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.
______________________________________________
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.
