R help - Jan 2012 - Function accepted by optim but not mle2 (?)

Dear Sir/ Madam,

I'm having trouble de-bugging the following - which works perfectly
well with optim or optimx - but not with mle2.
I'd be really grateful if someone could show me what is wrong. Many
thanks in advance. JSC:

gompertz<- function (x,t=data)
{
        a3<-x[1]
        b3<-x[2]
        shift<-data[1]
        h.t<-a3*exp(b3*(t-shift))
        S.t<-exp(a3/b3*(1-exp(b3*(t-shift))))
	PF<-1-(exp(a3/b3*(1-exp(b3*(data[1]-shift)))))
	QF<-1-(exp(a3/b3*(1-exp(b3*(data["length"(data)]-shift)))))
	return(S.t*h.t/(QF-PF))
}

## GOMPERTZ - maximum likelihood
gompertz2<-function (x)
{
lnlk<-sum(log(gompertz(x,t=data)))
return(-lnlk)
}

## this works:
A<-optimx(c(0.1,0.1),gompertz2,
method=c("Nelder-Mead"),control=list(fnscale=1))

## this doesn't:
A <- mle2(gompertz2(x),start=list(a3=0.03,b3=0.1),data=data)
## Error in gompertz2(x) : object 'x' not found

What is wrong is that you define and use parameters and global data
inconsistently in your functions.

v <- 6
f <- function( x ) { x^2 }
g <- function( x ) { v^3 }
f(v) # right answer
g(v) # right answer, accidentally
f(5) # works as expected
g(5) # surprise

You need to go through your functions and confirm for each variable that you
either define it in the function before you use it or that it is one of the
named parameters. It is possible to use global variables inside your functions,
but that is likely to give confusing answers so you should avoid it in most
cases (at least until you understand it).
---------------------------------------------------------------------------
Jeff Newmiller                        The     .....       .....  Go Live...
DCN:<jdnewmil at dcn.davis.ca.us>        Basics: ##.#.       ##.#.  Live
Go...
                                      Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
/Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k
--------------------------------------------------------------------------- 
Sent from my phone. Please excuse my brevity.

Medimel <medimelenglish at gmail.com> wrote:
>Dear Sir/ Madam,
>
>I'm having trouble de-bugging the following - which works perfectly
>well with optim or optimx - but not with mle2.
>I'd be really grateful if someone could show me what is wrong. Many
>thanks in advance. JSC:
>
>gompertz<- function (x,t=data)
>{
>        a3<-x[1]
>        b3<-x[2]
>        shift<-data[1]
>        h.t<-a3*exp(b3*(t-shift))
>        S.t<-exp(a3/b3*(1-exp(b3*(t-shift))))
>	PF<-1-(exp(a3/b3*(1-exp(b3*(data[1]-shift)))))
>	QF<-1-(exp(a3/b3*(1-exp(b3*(data["length"(data)]-shift)))))
>	return(S.t*h.t/(QF-PF))
>}
>
>## GOMPERTZ - maximum likelihood
>gompertz2<-function (x)
>{
>lnlk<-sum(log(gompertz(x,t=data)))
>return(-lnlk)
>}
>
>## this works:
>A<-optimx(c(0.1,0.1),gompertz2,
>method=c("Nelder-Mead"),control=list(fnscale=1))
>
>## this doesn't:
>A <- mle2(gompertz2(x),start=list(a3=0.03,b3=0.1),data=data)
>## Error in gompertz2(x) : object 'x' not found
>
>______________________________________________
>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.

Medimel <medimelenglish <at> gmail.com> writes:
> 
> Dear Sir/ Madam,
> 
> I'm having trouble de-bugging the following - which works perfectly
> well with optim or optimx - but not with mle2.
> I'd be really grateful if someone could show me what is wrong. Many
> thanks in advance. JSC:
> 


   mle2 is built on mle from the stats3 package, which expects
the function to be specified in terms of a set of independent parameters
(e.g.  gompertz2 <- function(a3,b3) {  ... }  

  There is a workaround, which could be better described in the documentation,
which is to set the 'vecpar' argument to TRUE and use parnames() a
parameter-names attribute on the function, as follows:

library(bbmle)

parnames(gompertz2) <- c("a3","b3")
A <- mle2(gompertz2,
          vecpar=TRUE,
          start=list(a3=0.03,b3=0.1),data=list(data=data))

Also note that 'data' must be specified as a *named* list,
as above.

  Your example was quite clear but not entirely reproducible --
it would have been good to provide a short example 'data' set.