Re: [R] [optim/bbmle] function returns NA

2013-08-14 Thread Ravi Varadhan
Carlos,
There are likely several problems with your likelihood.  You should check it 
carefully first before you do any optimization.  It seems to me that you have 
box constraints on the parameters.  They way you are enforcing them is not 
correct. I would prefer to use an optimization algorithm that can handle box 
constraints rather than use artificial mechanisms such as what you are trying 
to do:
if (r=0 | alpha=0 | s=0 | beta=0) return (NaN)

Even here, a better approach would be:
if (r=0 | alpha=0 | s=0 | beta=0) return (-.Machine$double.xmax)

I also notice that you do not use `tx' in your `g' function.  There are likely 
a number of other issues as well.

You may try the nmkb() function in the dfoptim package.  It can handle box 
constraints and typically tends to perform a bit better than optim's 
Nelder-Mead for unconstrained problems.

Hope this is helpful,
Ravi

[[alternative HTML version deleted]]

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


Re: [R] [optim/bbmle] function returns NA at

2013-08-13 Thread Prof J C Nash (U30A)
1) Why use Nelder-Mead with optimx when it is an optim() function. You 
are going from New York to Philadelphia via Beijing because of the extra 
overhead. The NM method is there for convenience in comparisons.


2) NM cannot work with NA when it wants to compute the centroid of 
points and search direction. So you've got to find a way to make sure 
your likelihood is properly defined. This seems to be the issue for 
about 90% of failures with optim(x) or other ML methods in my recent 
experience. Note that returning a large value (and make it a good deal 
smaller than the .Machine$double.xmax, say that number *1e-6 to avoid 
computation troubles) often works, but it is a quick and dirty fix.


JN


On 13-08-13 06:00 AM, r-help-requ...@r-project.org wrote:

Message: 36
Date: Tue, 13 Aug 2013 10:38:05 +0200
From: Carlos Nashercarlos.nas...@googlemail.com
To:r-help@r-project.org
Subject: [R] [optim/bbmle] function returns NA at ... distance from x
Message-ID:
CAP=bvwpxj991fbyt9ou5x1jf9nol3vtq1svtjvw82jwfjyz...@mail.gmail.com
Content-Type: text/plain

Dear R helpers,

I try to find the model parameters using mle2 (bbmle package). As I try to
optimize the likelihood function the following error message occurs:

Error in grad.default(objectivefunction, coef) :
   function returns NA at
1e-040.001013016911639890.0003166929388711890.000935163594829395 distance
from x.
In addition: Warning message:
In optimx(par = c(0.5, 10, 0.7, 10), fn = function (p)  :
   Gradient not computable after method Nelder-Mead

I can't figure out what that means exactly and how to fix it. I understand
that mle2 uses optim (or in my case optimx) to optimize the likelihood
function. As I use the Nelder-Mead method it should not be a problem if the
function returns NA at any iteration (as long as the initial values don't
return NA). Can anyone help me with that?

Here a small example of my code that reproduces the problem:

library(plyr)
library(optimx)

### Sample data ###
x - c(1,1,4,2,3,0,1,6,0,0)
tx - c(30.14, 5.14, 24.43, 10.57, 25.71, 0.00, 14.14, 32.86, 0.00, 0.00)
T - c(32.57, 29.14, 33.57, 34.71, 27.71, 38.14, 36.57, 37.71, 35.86, 30.57)
data - data.frame(x=x, tx=tx, T=T)

### Likelihood function ###
Likelihood - function(data, r, alpha, s, beta) {
   with(data, {
 if (r=0 | alpha=0 | s=0 | beta=0) return (NaN)
 f - function(x, tx, T)
 {
   g - function(y)
 (y + alpha)^(-( r + x))*(y + beta)^(-(s + 1))
   integrate(g, tx, T)$value
 }
 integral - mdply(data, f)
 L -
exp(lgamma(r+x)-lgamma(r)+r*(log(alpha)-log(alpha+T))-x*log(alpha+T)+s*(log(beta)-log(beta+T)))+exp(lgamma(r+x)-lgamma(r)+r*log(alpha)+log(s)+s*log(beta)+log(integral$V1))
 f - sum(log(L))
 return (f)
   })
}

### ML estimation function ###
Estimate_parameters_MLE - function(data, initValues) {
   llhd - function(r, alpha, s, beta) {
 return (Likelihood(data, r, alpha, s, beta))
   }
   library(bbmle)
   fit - mle2(llhd, initValues, skip.hessian=TRUE, optimizer=optimx,
method=Nelder-Mead, control=list(maxit=1e8))
   return (fit)
}

### Parameter estimation ###
Likelihood(data=data, r=0.5, alpha=10, s=0.7, beta=10) ### check initial
parameters -- -72.75183 -- initial parameters do return value
MLE_estimation - Estimate_parameters_MLE(data=data, list(r=0.5, alpha=10,
s=0.7, beta=10))

'Error in grad.default(objectivefunction, coef) :
   function returns NA at
1e-040.001013016911639890.0003166929388711890.000935163594829395 distance
from x.
In addition: Warning message:
   In optimx(par = c(0.5, 10, 0.7, 10), fn = function (p)  :
   Gradient not computable after method Nelder-Mead'


Best regards,
Carlos

-
Carlos Nasher
Buchenstr. 12
22299 Hamburg

tel:+49 (0)40 67952962
mobil:+49 (0)175 9386725
mail:carlos.nas...@gmail.com

[[alternative HTML version deleted]]



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