Christophe Pouzat
2005-Jul-21 16:01 UTC
[R] About object of class mle returned by user defined functions
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
Christophe Pouzat
2005-Jul-22 13:09 UTC
[R] About object of class mle returned by user defined functions
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