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