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 <[EMAIL PROTECTED]>
To: r-help@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.
______________________________________________
R-help@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.
