[R] About object of class mle returned by user defined functions

Christophe Pouzat christophe.pouzat at univ-paris5.fr
Fri Jul 22 15:09:30 CEST 2005


Guys,

I apologize for being slightly misleading in my previous e-mail.

First, I generated some confusion between the scale and rate parameters 
in the gamma distribution. My direct call to mle use a minuslogl 
function "working" with a scale parameter while my call to mle from my 
function used a minuslogl function "working" with a rate parameter!... 
To add to the confusion I had simulated data with a scale ( = 1/rate) 
value of 1... I really hope that none of you lost time with that.

Second, some "^2" in my original function definition got converted into 
exponents on the e-mail, meaning that if some of you tried to copy and 
paste it you must have gotten some insults from R while sourcing it. In 
order to avoid that I attach an ".R" file. In principle if you source it 
and then type the following commands you should get (exactly):


 > coef(fitA)
    shape     scale
2.2230421 0.8312374
 > coef(fit1)
    shape     scale
2.2230421 0.8312374
 > vcov(fitA)
            shape       scale
shape  0.08635158 -0.03228829
scale -0.03228829  0.01518126
 > vcov(fit1)
            shape       scale
shape  0.08635158 -0.03228829
scale -0.03228829  0.01518126
 > logLik(fitA)
'log Lik.' -146.6104 (df=2)
 > logLik(fit1)
'log Lik.' -146.6104 (df=2)
 > confint(fitA)
Profiling...
          2.5 %   97.5 %
shape 1.6985621 2.853007
scale 0.6307824 1.129889
 > confint(fit1)
Profiling...
Erreur dans approx(sp$y, sp$x, xout = cutoff) :
    need at least two non-NA values to interpolate
De plus : Message d'avis :
collapsing to unique 'x' values in: approx(sp$y, sp$x, xout = cutoff)

Here fitA is obtained by a direct call to mle (I mean from the command 
line) while fit1 is obtained by the same call but within a function: newFit.

The fundamental problem remains, I don't understand why confint does 
work with fitA and not with fit1.

Christophe.

PS: my version info

platform i686-pc-linux-gnu
arch     i686            
os       linux-gnu       
system   i686, linux-gnu 
status                   
major    2               
minor    1.1             
year     2005            
month    06              
day      20              
language R



Christophe Pouzat wrote:

>Hi,
>
>There is something I don't get with object of class "mle" returned by a 
>function I wrote. More precisely it's about the behaviour of method 
>"confint" and "profile" applied to these object.
>
>I've written a short function (see below) whose arguments are:
>1) A univariate sample (arising from a gamma, log-normal or whatever).
>2) A character string standing for one of the R densities, eg, "gamma", 
>"lnorm", etc. That's the density the user wants to fit to the data.
>3) A named list with initial values for the density parameters; that 
>will be passed to optim via mle.
>4) The method to be used by optim via mle. That can be change by the 
>code if parameter boundaries are also supplied.
>5) The lowest allowed values for the parameters.
>6) The largest allowed values.
>
>The "big" thing this short function does is writing on-fly the 
>corresponding log-likelihood function before calling "mle". The object 
>of class "mle" returned by the call to "mle" is itself returned by the 
>function.
>
>Here is the code:
>
>newFit <- function(isi, ## The data set
>                   isi.density = "gamma", ## The name of the density 
>used as model
>                   initial.para = list( shape = (mean(isi)/sd(isi))^2,
>                     scale = sd(isi)^2 / mean(isi) ), ## Inital 
>parameters passed to optim
>                   optim.method = "BFGS", ## optim method
>                   optim.lower = numeric(length(initial.para)) + 0.00001,
>                   optim.upper = numeric(length(initial.para)) + Inf,
>                   ...) {
>
>  require(stats4)
> 
>  ## Create a string with the log likelihood definition
>  minusLogLikelihood.txt <- paste("function( ",
>                                  paste(names(initial.para), collapse = 
>", "),
>                                  " ) {",
>                                  "isi <- eval(",
>                                  deparse(substitute(isi)),
>                                  ", envir = .GlobalEnv);",
>                                  "-sum(",
>                                  paste("d", isi.density, sep = ""),
>                                  "(isi, ",
>                                  paste(names(initial.para), collapse = 
>", "),
>                                  ", log = TRUE) ) }"
>                                  )
>
>  ## Define logLikelihood function
>  minusLogLikelihood <- eval( parse(text = minusLogLikelihood.txt) )
>  environment(minusLogLikelihood) <- .GlobalEnv
>
> 
>  if ( all( is.infinite( c(optim.lower,optim.upper) ) ) ) {
>      getFit <- mle(minusLogLikelihood,
>                    start = initial.para,
>                    method = optim.method,
>                    ...
>                    )
>  } else {
>    getFit <- mle(minusLogLikelihood,
>                  start = initial.para,
>                  method = "L-BFGS-B",
>                  lower = optim.lower,
>                  upper = optim.upper,
>                  ...
>                  )
>  }  ## End of conditional on all(is.infinite(c(optim.lower,optim.upper)))
> 
>  getFit
> 
>}
>
>
>It seems to work fine on examples like:
>
> > isi1 <- rgamma(100, shape = 2, scale = 1)
> > fit1 <- newFit(isi1) ## fitting here with the "correct" density 
>(initial parameters are obtained by the method of moments)
> > coef(fit1)
>    shape     scale
>1.8210477 0.9514774
> > vcov(fit1)
>           shape      scale
>shape 0.05650600 0.02952371
>scale 0.02952371 0.02039714
> > logLik(fit1)
>'log Lik.' -155.9232 (df=2)
>
>If we compare with a "direct" call to "mle":
>
> > llgamma <- function(sh, sc) -sum(dgamma(isi1, shape = sh, scale = sc, 
>log = TRUE))
> > fitA <- mle(llgamma, start = list( sh = (mean(isi1)/sd(isi1))^2, sc = 
>sd(isi1)^2 / mean(isi1) ),lower = c(0.0001,0.0001), method = "L-BFGS-B")
> > coef(fitA)
>      sh       sc
>1.821042 1.051001
> > vcov(fitA)
>            sh          sc
>sh  0.05650526 -0.03261146
>sc -0.03261146  0.02488714
> > logLik(fitA)
>'log Lik.' -155.9232 (df=2)
>
>I get almost the same estimated parameter values, same log-likelihood 
>but not the same vcov matrix.
>
>A call to "profile" or "confint" on fit1 does not work, eg:
> > confint(fit1)
>Profiling...
>Erreur dans approx(sp$y, sp$x, xout = cutoff) :
>    need at least two non-NA values to interpolate
>De plus : Message d'avis :
>collapsing to unique 'x' values in: approx(sp$y, sp$x, xout = cutoff)
>
>Although calling the log-likelihood function defined in fit1 
>(fit1 at minuslogl) with argument values different from the MLE does return 
>something sensible:
>
> > fit1 at minuslogl(coef(fit1)[1],coef(fit1)[2])
>[1] 155.9232
> > fit1 at minuslogl(coef(fit1)[1]+0.01,coef(fit1)[2]+0.01)
>[1] 155.9263
>
>There is obviously something I'm missing here since I thought for a 
>while that the problem was with the environment "attached" to the 
>function "minusLogLikelihood" when calling "eval"; but the lines above 
>make me think it is not the case...
>
>Any help and/or ideas warmly welcomed.
>
>Thanks,
>
>Christophe.
>
>  
>


-- 
A Master Carpenter has many tools and is expert with most of them.If you
only know how to use a hammer, every problem starts to look like a nail.
Stay away from that trap.
Richard B Johnson.
--

Christophe Pouzat
Laboratoire de Physiologie Cerebrale
CNRS UMR 8118
UFR biomedicale de l'Universite Paris V
45, rue des Saints Peres
75006 PARIS
France

tel: +33 (0)1 42 86 38 28
fax: +33 (0)1 42 86 38 30
web: www.biomedicale.univ-paris5.fr/physcerv/C_Pouzat.html

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: testScript.R
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20050722/50371b4c/testScript.pl


More information about the R-help mailing list