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]]