[R] Likelihood optimization numerically
Bill.Venables at csiro.au
Bill.Venables at csiro.au
Sun Jan 27 08:13:18 CET 2008
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:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of Mohammad Ehsanul Karim
Sent: Sunday, 27 January 2008 4:48 PM
To: r-help at 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 at 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.
More information about the R-help
mailing list