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@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.