[R] Likelihood optimization numerically

[Bill.Venables at csiro.au]

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.