[R] Error in optim when i call it from a function

From: Borja Soto Varela <borja.soto_at_gmail.com>
Date: Fri, 14 Nov 2008 21:12:55 +0100


Dear R-users

I've got the next problem:

I've got this *function*:

fitcond=function(x,densfun,pcorte,start,...){

    myfn <- function(parm,x,pcorte,...) -sum(log(dens(parm,x,pcorte,...)))     Call <- match.call(expand.dots = TRUE)     if (missing(start))

        start <- NULL
    dots <- names(list(...))
    dots <- dots[!is.element(dots, c("upper", "lower"))]     if (missing(x) || length(x) == 0 || mode(x) != "numeric")

        stop("'x' must be a non-empty numeric vector")     if (missing(densfun) || !(is.function(densfun) || is.character(densfun)))

        stop("'densfun' must be supplied as a function or name")     n <- length(x)
    if (is.character(densfun)) {

        distname <- tolower(densfun)
        densfun <- switch(distname,exponential = dexp,gamma = dgamma,
"log-normal" = dlnorm,
            lognormal = dlnorm, weibull = dweibull, pareto = dpareto,
            NULL)
      distfun <- switch(distname,exponential = pexp,gamma = pgamma,
"log-normal" = plnorm,
            lognormal = plnorm, weibull = pweibull, pareto = ppareto,
            NULL)

        if (is.null(densfun))
            stop("unsupported distribution")
        if (distname %in% c("lognormal", "log-normal")) {
            if (!is.null(start))
                stop("supplying pars for the log-Normal is not supported")
            if (any(x <= 0))
                stop("need positive values to fit a log-Normal")
            lx <- log(x)
            sd0 <- sqrt((n - 1)/n) * sd(lx)
            mx <- mean(lx)
        start <- list(meanlog = mx, sdlog = sd0)
            start <- start[!is.element(names(start), dots)]
        }
        if (distname == "exponential") {
            if (!is.null(start))
                stop("supplying pars for the exponential is not supported")
            estimate <- 1/mean(x)
            #sds <- estimate/sqrt(n)
            names(estimate) <- names(sds) <- "rate"
        start <- list(rate=estimate)
        start <- start[!is.element(names(start),dots)]
        }
        if (distname == "weibull" && is.null(start)) {
            m <- mean(log(x))
            v <- var(log(x))
            shape <- 1.2/sqrt(v)
            scale <- exp(m + 0.572/shape)
            start <- list(shape = shape, scale = scale)
            start <- start[!is.element(names(start), dots)]
        }
        if (distname == "gamma" && is.null(start)) {
            m <- mean(x)
            v <- var(x)
            start <- list(shape = m^2/v, rate = m/v)
            start <- start[!is.element(names(start), dots)]
        }

    }
    if (is.null(start) || !is.list(start))

        stop("'start' must be a named list")     nm <- names(start)
    f <- formals(densfun)
    args <- names(f)
    m <- match(nm, args)
    if (any(is.na(m)))

        stop("'start' specifies names which are not arguments to 'densfun'")     formals(densfun) <- c(f[c(1, m)], f[-c(1, m)])     #dens <- function(parm, x,pcorte, ...) densfun(x, parm,
...)/(1-distfun(pcorte,parm))

    #if ((l <- length(nm)) > 1)
    #   body(dens) <- parse(text = paste("densfun(x,", paste("parm[",
    #        1:l, "]", collapse = ", "), ",
...)","/(1-distfun(pcorte,",paste("parm[",
    #       1:l, "]", collapse = ", "), ", ...))"))

    Call[[1]] <- as.name("optim")

    Call$densfun <- Call$start <- NULL
    Call$x <- x
    Call$pcorte <- pcorte
    Call$par <- start
    Call$hessian <- TRUE

    fn1=function(parm,x,pcorte,...) -sum(log(dgamma(x, parm[1], parm[2],
...)/(1 - pgamma(pcorte, parm[1],parm[2], ... ))))

    Call$fn <- fn1
    if (is.null(Call$method)) {

       if (any(c("lower", "upper") %in% names(Call)))
            Call$method <- "L-BFGS-B"
        else if (length(start) > 1)
           Call$method <- "BFGS"
        else Call$method <- "Nelder-Mead"
    }
    print(Call)
    res <- eval.parent(Call)
    list(estimate = res$par, loglik = -res$value) }

*When I call it , I've got the next problem:* >corte=10.51504

>datos1=c(10.93822, 11.74718, 10.86152, 11.14968, 11.34962,
11.00662, 10.98304, 11.06032, 11.33369, 12.61333, 10.55254, 10.98838,
10.61359, 10.54633, 11.48311, 12.24803, 10.98632, 10.53006, 12.13353,
10.93739, 10.85004, 10.95945, 13.70617, 10.785, 11.2085, 11.8517,
10.99292, 11.01403, 10.59838, 13.92816, 10.81108, 11.14236, 10.82758,
13.07142)
>fitcond(datos1,"gamma",pcorte=10.51504,control=list(maxit=1000))

optim(x = c(10.93822, 11.74718, 10.86152, 11.14968, 11.34962,

11.00662, 10.98304, 11.06032, 11.33369, 12.61333, 10.55254, 10.98838,
10.61359, 10.54633, 11.48311, 12.24803, 10.98632, 10.53006, 12.13353,
10.93739, 10.85004, 10.95945, 13.70617, 10.785, 11.2085, 11.8517,
10.99292, 11.01403, 10.59838, 13.92816, 10.81108, 11.14236, 10.82758,
13.07142), pcorte = 10.51504, control = list(maxit = 1000), par = list(

    shape = 173.623466051441, rate = 15.3008183026100), hessian = TRUE,     fn = function(parm,x,pcorte,...) -sum(log(dgamma(x, parm[1], parm[2],
...)/(1 - pgamma(pcorte, parm[1],parm[2], ... ))))

    , method = "BFGS")
Error en optim(x = c(10.93822, 11.74718, 10.86152, 11.14968, 11.34962, :   non-finite finite-difference value [1] Además: Hubo 30 avisos (use warnings() para verlos)

but if i try only the call to optim like I print inside the function it runs properly:

optim(x = c(10.93822, 11.74718, 10.86152, 11.14968, 11.34962,

11.00662, 10.98304, 11.06032, 11.33369, 12.61333, 10.55254, 10.98838,
10.61359, 10.54633, 11.48311, 12.24803, 10.98632, 10.53006, 12.13353,
10.93739, 10.85004, 10.95945, 13.70617, 10.785, 11.2085, 11.8517,
10.99292, 11.01403, 10.59838, 13.92816, 10.81108, 11.14236, 10.82758,
13.07142), pcorte = 10.51504, control = list(maxit = 1000), par = list(

    shape = 173.623466051441, rate = 15.3008183026100), hessian = TRUE,     fn = function(parm,x,pcorte,...) -sum(log(dgamma(x, parm[1], parm[2],
...)/(1 - pgamma(pcorte, parm[1],parm[2], ... ))))

    , method = "BFGS").
$par

    shape rate
0.0035234 1.1185760

Does someone know what is happening?

Thanks
Borja

        [[alternative HTML version deleted]]



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Fri 14 Nov 2008 - 20:15:38 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 14 Nov 2008 - 21:30:25 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive