[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