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