[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) :
onediml optimization by NelderMead 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: rhelp at rproject.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
loglikelihood
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]
"i386pcmingw32"
$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
(20071126)"
____________________________________________________________________________________
Never miss a thing. Make Yahoo your home page.
More information about the Rhelp
mailing list