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

From: Christophe Pouzat <christophe.pouzat_at_univ-paris5.fr>
Date: Fri 22 Jul 2005 - 23:09:30 EST

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@minuslogl) with argument values different from the MLE does return
>something sensible:
>
> > fit1@minuslogl(coef(fit1)[1],coef(fit1)[2])
>[1] 155.9232
> > fit1@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


set.seed(121, "Mersenne-Twister") ## Set the rng seed ## We will need stats4 library library(stats4) ## Define function newFit 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) para.names <- paste(names(initial.para)[1], " = ", names(initial.para)[1], sep = "") if (length(names(initial.para)) > 1) for ( this.name in names(initial.para)[-1] ) { para.names <- paste(para.names, ", ", this.name, " = ", this.name, sep = "" ) } ## Create a string with the log likelihood definition minusLogLikelihood.txt <- paste("function( ", paste(names(initial.para), collapse = ", "), " )", "-sum(", paste("d", isi.density, "(", deparse(substitute(isi)), ",", para.names, ", log = TRUE) )", sep = "" ) ) ## 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 } ## End of function newFit definition isi1 <- rgamma(100, shape = 2, scale = 1) ## Generate 100 rn with a Gamma distribution ## Define function llgama which return the oposite of the log-likelihood of isi1 assuming ## a Gamma density with parameters shape = shape and scale = scale llgamma <- function(shape, scale) -sum(dgamma(isi1, shape = shape, scale = scale, log = TRUE)) ## Call mle to perform the fit using the method of moments for the initial guess fitA <- mle(llgamma, start = list( shape = (mean(isi1)/sd(isi1))^2, scale = sd(isi1)^2 / mean(isi1) ),lower = c(0.00001,0.00001), method = "L-BFGS-B") ## Do in principle the same thing by calling newFit fit1 <- newFit(isi = isi1, isi.density = "gamma", initial.para = list( shape = (mean(isi1)/sd(isi1))^2, scale = sd(isi1)^2 / mean(isi1) ), optim.method = "L-BFGS-B", optim.lower = c(0.00001,0.00001))

______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Received on Fri Jul 22 23:17:11 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:33:55 EST