# [R] Likelihood optimization numerically

Mohammad Ehsanul Karim wildscop at yahoo.com
Sun Jan 27 10:11:18 CET 2008

```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
time.

Ehsan

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

____________________________________________________________________________________