Could you use 'optim' directly? This won't return a nice object of
class 'mle', but it will give you parameter estimates.
If you want an object of class 'mle', have you worked through the
examples in the help file?
Also, are your arguments "a" and "b" are scalars? If yes, did you
try passing them to your 'minuslogl' function in the 'fixed' argument?
If this won't work for you, you can try the functions 'mle.' and
'profile.mle' below. If I'm not mistaken, the standard 'mle' requires
the 'fixed' argument to be a named list of scalars. I don't remember
now for sure, but I think the 'mle.' function below differs from the
standard 'mle' only in allowing 'fixed' arguments of length greater than
1 and then dropping them from operations that require scalars. My
'profile.mle' is modified from the function in 'stats4' to make error
termination less likely, at least for the applications I used.
Hope this helps.
Spencer Graves
p.s. I thought about sending this to someone, but the "Maintainer" for
the "stats4" package is the "R Core Team <[EMAIL PROTECTED]>", and I
decided not to bother that whole group.
#######################################
mle. <-
function(minuslogl, start = formals(minuslogl), method = "BFGS",
fixed = list(), ...)
{
call <- match.call()
n <- names(fixed)
fullcoef <- formals(minuslogl)
if (any(!n %in% names(fullcoef)))stop(
"some named arguments in 'fixed' are not arguments to the
supplied log-likelihood")
fullcoef[n] <- fixed
if (!missing(start) && (!is.list(start) || is.null(names(start))))
stop("'start' must be a named list")
start[n] <- NULL
start <- sapply(start, eval.parent)
nm <- names(start)
oo <- match(nm, names(fullcoef))
if (any(is.na(oo)))
stop("some named arguments in 'start' are not arguments to the
supplied log-likelihood")
start <- start[order(oo)]
f <- function(p) {
l <- as.list(p)
names(l) <- nm
l[n] <- fixed
do.call("minuslogl", l)
}
dots <- list(...)
oout <- optim(start, f, method = method, hessian = TRUE,
...)
oout$fixed <- fixed
oout$dots <- list(...)
coef <- oout$par
vcov <- if (length(coef))
solve(oout$hessian)
else matrix(numeric(0), 0, 0)
min <- oout$value
fcoef <- fullcoef
fcoef[nm] <- coef
len.coef <- sapply(fcoef, length)
# Drop "fixed" arguments of length != 1 as they
# almost certainly relate arguments of "minuslogl"
# that are NOT parameters to be estimated.
f.coef <- fcoef[len.coef==1]
newMLE <- new("mle", call = call, coef = coef, fullcoef =
unlist(f.coef),
vcov = vcov, min = min, details = oout, minuslogl = minuslogl,
method = method)
}
#########################
#getMethods("profile", "package:stats4")
profile.mle <- function (fitted, ...){
.local <- function (fitted, which = 1:p, maxsteps = 100,
alpha = 0.01, zmax = sqrt(qchisq(1 - alpha/2, p)), del = zmax/5,
trace = FALSE, ...)
{
onestep <- function(step, call, fixed) {
bi <- B0[i] + sgn * step * del * std.err[i]
# Additional fix for profile
# Get original fix
{
if(is.list(fixed) && (length(fixed)>0)){
fixed[[pi]] <- NULL
k. <- length(fixed)
fix2 <- vector(k.+1, mode="list")
names(fix2) <- c(names(fixed), pi)
fix2[1:k.] <- fixed
fix2[k.+1] <- bi
}
else{
fix2 <- list(bi)
names(fix2) <- pi
}
}
call$fixed <- fix2
pfit <- try(eval.parent(call, 2), silent = TRUE)
if (inherits(pfit, "try-error"))
return(NA)
else {
zz <- 2 * ([EMAIL PROTECTED] - [EMAIL PROTECTED])
ri <- pv0
ri[, names([EMAIL PROTECTED])] <- [EMAIL PROTECTED]
ri[, pi] <- bi
if(zz<=0){
msg0 <- paste("Profiling has found a better ",
"solution (log(LR) = ", zz, " <0), so the ",
"original fit had not converged: ",
"This will be ignored for now, but try ",
"mle(..., start=list(", sep="")
pars <- paste(Pnames, coef(pfit), sep="=")
pars. <- paste(pars, collapse=", ")
msg1 <- paste(msg0, pars., ")).\n", sep="")
warning(msg1)
}
z0 <- max(zz, 0)
z <- sgn * sqrt(z0)
pvi <<- rbind(pvi, ri)
zi <<- c(zi, z)
}
if (trace)
cat(bi, z, "\n")
########### end subfunction onestep ##########
list(z=z, zfit=pfit)
}
betterSum <- function(fit){
if((!inherits(fit, "try-error")) && (class(fit)=="mle")){
fitMin <- [EMAIL PROTECTED]
if(!is.numeric(fitMin)){
warn("'mle' returned 'min' of class ", class(fitMin),
"; must be numeric.")
return(B1.)
}
if((nMin <- length(fitMin))!=1){
warn("'mle' returned 'min' of length ", nMin,
"; must be 1")
return(B1.)
}
if(fitMin<B1.$minuslogl)
return(list([EMAIL PROTECTED],[EMAIL PROTECTED],
nImprovementsFound=B1.$nImprovementsFound+1))
}
########### end subfunction betterSum ##########
B1.
}
summ <- summary(fitted)
std.err <- [EMAIL PROTECTED], "Std. Error"]
Pnames <- names(B0 <- [EMAIL PROTECTED])
pv0 <- t(as.matrix(B0))
B0. <- list(coef=B0, [EMAIL PROTECTED])
B1. <- list(coef=B0, [EMAIL PROTECTED],
nImprovementsFound=0)
p <- length(Pnames)
prof <- vector("list", length = length(which))
names(prof) <- Pnames[which]
call <- [EMAIL PROTECTED]
call$minuslogl <- [EMAIL PROTECTED]
call$method <- [EMAIL PROTECTED]
ndeps <- eval.parent(call$control$ndeps)
parscale <- eval.parent(call$control$parscale)
fixed <- [EMAIL PROTECTED]
dots <- [EMAIL PROTECTED]
if(length(dots)>0)
for(i.. in names(dots))
call[[i..]] <- dots[[i..]]
for (i in which) {
zi <- 0
pvi <- pv0
pi <- Pnames[i]
if (!is.null(ndeps))
call$control$ndeps <- ndeps[-i]
if (!is.null(parscale))
call$control$parscale <- parscale[-i]
for (sgn in c(-1, 1)) {
if (trace)
cat("\nParameter:", pi, c("down", "up")[(sgn +
1)/2 + 1], "\n")
step <- 0
z <- 0
call$start <- as.list(B0)
lastz <- 0
while ((step <- step + 1) < maxsteps && abs(z) <
zmax) {
z. <- onestep(step, call, fixed)
B1. <- betterSum(z.$zfit)
z <- z.$z
if (is.na(z))
break
lastz <- z
}
if (abs(lastz) < zmax) {
for (dstep in c(0.2, 0.4, 0.6, 0.8, 0.9)) {
z. <- onestep(step-1+dstep, call, fixed)
B1. <- betterSum(z.$zfit)
z <- z.$z
if (is.na(z) || abs(z) > zmax)
break
}
}
else if (length(zi) < 5) {
mxstep <- step - 1
step <- 0.5
while ((step <- step + 1) < mxstep){
z. <- onestep(step, call, fixed)
B1. <- betterSum(z.$zfit)
# z <- z.$z
}
}
}
si <- order(pvi[, i])
prof[[pi]] <- data.frame(z = zi[si])
prof[[pi]]$par.vals <- pvi[si, , drop = FALSE]
}
new("profile.mle", profile = prof, summary = summ)
}
.local(fitted, ...)
}
###################################################
# mle examples:
x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
ll <- function(ymax=15, xhalf=6)
-sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE))
(fit <- mle(ll))
(fit. <- mle.(ll))
(fit.5 <- mle.(ll, fixed=list(xhalf=6)))
summary(fit)
summary(fit.)
logLik(fit)
logLik(fit.)
vcov(fit)
vcov(fit.)
plot(profile(fit), absVal=FALSE)
plot(prof. <- profile.mle(fit.), absVal=FALSE)
confint(fit)
confint(prof.)
## use bounded optimization
## the lower bounds are really > 0, but we use >=0 to stress-test
profiling
(fit1 <- mle(ll, method="L-BFGS-B", lower=c(0, 0)))
(fit1. <- mle.(ll, method="L-BFGS-B", lower=c(0, 0)))
plot(profile(fit1), absVal=FALSE)
plot(profile.mle(fit1), absVal=FALSE)
## a better parametrization:
ll2 <- function(lymax=log(15), lxhalf=log(6))
-sum(stats::dpois(y, lambda=exp(lymax)/(1+x/exp(lxhalf)),
log=TRUE))
(fit2 <- mle(ll2))
(fit2. <- mle.(ll2))
plot(profile(fit2), absVal=FALSE)
plot(prof2 <- profile.mle(fit2), absVal=FALSE)
exp(confint(fit2))
exp(confint(prof2))
ll. <- function(ymax=15, xhalf=6, x, y)
-sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE))
fit <- mle(ll)
fit. <- mle.(ll)
all.equal(fit, fit.)
prof <- profile(fit)
plot(prof)
with([EMAIL PROTECTED],
plot(par.vals[,2], z, type="l"))
with([EMAIL PROTECTED],
plot(par.vals[,2], abs(z), type="l"))
with([EMAIL PROTECTED],
plot(par.vals[,1], abs(z), type="l"))
undebug(profile.mle)
prof. <- profile.mle(fit.)
plot(prof.)
(fit.xy <- mle.(ll., start=list(ymax=20, xhalf=3),
fixed=list(x=x, y=y)))
(prof.fit <- profile.mle(fit.xy))
plot(prof.fit)
logLik(fit.xy)
logLik(fit)
logLik(fit.)
mle2 <- function(x, y){
mle.(ll., start=list(ymax=20, xhalf=3),
fixed=list(x=x, y=y))
}
fit2 <- mle2(x, y)
prof2 <- profile.mle(fit2)
plot(prof2)
confint(prof2)
mle3 <- function(x, y){
do.call('mle.',
list(ll., start=list(ymax=20, xhalf=3),
fixed=list(x=x, y=y)))
}
fit3 <- mle3(x, y)
prof3 <- profile.mle(fit3)
confint(prof3)
plot(prof3)
mle4 <- function(x, y){
eval(mle.(ll., start=list(ymax=20, xhalf=3),
fixed=list(x=x, y=y)))
}
fit4 <- mle4(x, y)
prof4 <- profile.mle(fit4)
confint(prof4)
mle5 <- function(x, y){
mle.call <- paste(
"mle.(ll.,",
"start=list(ymax=20, xhalf=3),",
"fixed=list(x=x, y=y))")
mle.parsed <- parse(text=mle.call)
eval(mle.parsed)
}
debug(mle5)
fit5 <- mle5(x, y)
prof5 <- profile.mle(fit5)
confint(prof5)
plot(prof5)
#######################################
Rainer M Krug wrote:
> Hi
>
> I have the following formula (I hope it is clear - if no, I can try to
> do better the next time)
>
> h(x, a, b) =
> integral(0 to pi/2)
> (
> (
> integral(D/sin(alpha) to Inf)
> (
> (
> f(x, a, b)
> )
> dx
> )
> dalpha
> )
>
> and I want to do an mle with it.
> I know how to use mle() and I also know about integrate(). My problem is
> to give the parameter values a and b to the integrate function.
>
> In other words, how can I write
>
> h <- function...
>
> so that I can estimate a and b?
>
> Thanks,
>
> Rainer
>
>
______________________________________________
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.html
