[R] About object of class mle returned by user defined functions
Christophe Pouzat
christophe.pouzat at univ-paris5.fr
Thu Jul 21 18:01:36 CEST 2005
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
More information about the R-help
mailing list