[R] Question regarding to maxNR
Bill.Venables at csiro.au
Bill.Venables at csiro.au
Fri Mar 12 07:28:59 CET 2010
PS there is a better option still. Replace
log(prod(dcauchy(x,mu,s)))
with
sum(dcauchy(x,mu,s, log = TRUE))
For huge samples this will me milliseconds faster...
Bill Venables
CSIRO/CMIS Cleveland Laboratories
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Bill.Venables at csiro.au
Sent: Friday, 12 March 2010 4:23 PM
To: yhsu6 at illinois.edu; r-help at r-project.org
Subject: [ExternalEmail] Re: [R] Question regarding to maxNR
Your problem is numerical.
Try replacing
log(prod(dcauchy(x,mu,s)))
by
sum(log(dcauchy(x,mu,s)))
and see the difference. Here's what I get:
> mu <- 2
> s <- 1
> n <- 300
> library(maxLik)
> set.seed(1004)
> x <- rcauchy(n,mu,s)
> loglik <- function(mu) {
+ sum(log(dcauchy(x,mu,s)))
+ }
> maxNR(loglik,start=median(x))$estimate
[1] 2.075059
which looks about right to me.
Bill Venables
CSIRO/CMIS Cleveland Laboratories
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of kate
Sent: Friday, 12 March 2010 3:21 PM
To: r-help at r-project.org
Subject: [R] Question regarding to maxNR
Hi R-users,
Recently, I use maxNR function to find maximizer. I have error appears as follows
Error in maxNRCompute(fn = fn, grad = grad, hess = hess, start = start, :
NA in the initial gradient
My code is
mu=2
s=1
n=300
library(maxLik)
set.seed(1004)
x<-rcauchy(n,mu,s)
loglik<-function(mu)
{
log(prod(dcauchy(x,mu,s)))
}
maxNR(loglik,start=median(x))$estimate
Does anyone know how to solve this problem?
Thanks,
Kate
[[alternative HTML version deleted]]
______________________________________________
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.
______________________________________________
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.
More information about the R-help
mailing list