# 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,

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),coef(fit1))
> 155.9232
> > fit1@minuslogl(coef(fit1)+0.01,coef(fit1)+0.01)
> 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), " = ", names(initial.para), 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