[R] Likelihood optimization numerically

[Mohammad Ehsanul Karim]

Dear List,

Probably i am missing something important in optimize:

llk.1st <- function(alpha){
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)
llk1 <- -n*log(gamma(alpha)) - n*alpha*log(sum(x)/(n*alpha)) + (alpha - 1)*(sum(log(x))) - (sum(x))/(sum(x)/(n*alpha))
return(llk1)
} 

plot(llk.1st, 1,1000)
optimize(f=llk.1st, interval =c(0,1000), tol = 0.0001)


Reported result is approximately -
alpha = 500
beta = 0.05

But i get -
$minimum
[1] 1000 ## whatever the upper bound is!!
$objective
[1] -Inf
There were 50 or more warnings (use warnings() to see the first 50)

also,
> optim(par = 500, fn = llk.1st)
Error in optim(par = 500, fn = llk.1st) : 
  function cannot be evaluated at initial parameters
In addition: Warning messages:
1: In optim(par = 500, fn = llk.1st) :
  one-diml optimization by Nelder-Mead is unreliable: use optimize
2: In fn(par, ...) : value out of range in 'gammafn'


Can 
anyone 
provide 
me 
hint?
Thank 
you 
for 
your 
time.

Ehsan
http://www.youtube.com/profile_play_list?user=wildsc0p


----- Original Message ----
From: Mohammad Ehsanul Karim <wildscop at yahoo.com>
To: r-help at r-project.org
Sent: Sunday, January 27, 2008 12:47:40 AM
Subject: 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".

> 
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)"





      ____________________________________________________________________________________
Never miss a thing.  Make Yahoo your home page.