From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>

Date: Fri 22 Jul 2005 - 23:47:39 EST

*> fit1@call
*

mle(minuslogl = minusLogLikelihood, start = initial.para, method =**"L-BFGS-B",
**

lower = optim.lower, upper = optim.upper)

Date: Fri 22 Jul 2005 - 23:47:39 EST

*> fitA@call
*

mle(minuslogl = llgamma, start = list(shape = (mean(isi1)/sd(isi1))^2,

scale = sd(isi1)^2/mean(isi1)), method = "L-BFGS-B", lower = c(1e-05, 1e-05))

mle(minuslogl = minusLogLikelihood, start = initial.para, method =

lower = optim.lower, upper = optim.upper)

and the difference should be clear.

You can fix this up by constructing a call containing the values and not the names (one way is to use substitute()) and then eval() it. Another way is something like

Call <- quote(mle())

Call$minuslogl <- minusLogLikelihood

Call$start <- initial.para

...

On Fri, 22 Jul 2005, Christophe Pouzat wrote:

> 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
**>
**>
*

-- Brian D. Ripley, ripley@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595 ______________________________________________ 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.htmlReceived on Fri Jul 22 23:54:46 2005

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