diff --git a/R/fitmand.R b/R/fitmand.R index 931e730..0b4b011 100644 --- a/R/fitmand.R +++ b/R/fitmand.R @@ -8,7 +8,9 @@ fitmand <- function(x, trunc, start.value, ...){ if (min(y)<=trunc) stop("truncation point should be lower than the lowest data value") } if(missing(start.value)){ - shat <- 2 + x75 <- rad.tab[(floor(dim(rad.tab)[1]/4)):dim(rad.tab)[1],] + x75 <- x75[x75$abund > 1, ] + shat <- - coef(lm(log(abund)~log(rank), data=x75))[2] vhat <- 30 } else{ diff --git a/R/plotprofmle.R b/R/plotprofmle.R index 6807b2b..40530ce 100644 --- a/R/plotprofmle.R +++ b/R/plotprofmle.R @@ -41,6 +41,7 @@ plotprofmle <- function(profobj, nseg=20, ratio=log(8), which=1:length(profobj@p change <- (interpol$y - ratio)[2:l] * (interpol$y - ratio)[1:(l-1)] endpoints <- which(change < 0) corr <- (interpol$x[2]-interpol$x[1])/2 + if(length(endpoints) > 0) for (j in 1:(length(endpoints)/2)) { lower <-interpol$x[endpoints[(2*j)-1]]+corr upper <- interpol$x[endpoints[2*j]]+corr diff --git a/R/updatesad.R b/R/updatesad.R index eb607cb..1eaad8a 100644 --- a/R/updatesad.R +++ b/R/updatesad.R @@ -1,7 +1,7 @@ ## updates a fitsad/fitrad object by running the optimizer again starting ## on better fit returned by profile updatesad <- function(object, ...) { - dots <- as.list(...) + dots <- list(...) prof <- profile(object) if(class(prof) == "profile.mle2") stop("Cannot update, profile did not find a better fit!") newcall <- as.list(object@call) @@ -9,9 +9,15 @@ updatesad <- function(object, ...) { newcall[["control"]] <- NULL # removes the "control" slot # merges the original call with arguments supplied by ... for (v in names(dots)) newcall[[v]] <- dots[[v]] - names(newcall$start) -> name # bugfix? profile returns coefficient "prob" name as "prob.prob" - newcall$start <- list(prof@fullcoef) # TODO: might be coef instead? how to deal with fixed? - names(newcall$start) <- name + # bugfix? profile returns coefficient "prob" name as "prob.prob", sometimes changes parameters + # from start to fixed... + names(newcall$start) -> name + for (v in name) { + if (v %in% names(prof@coef)) + newcall$start[[v]] <- as.numeric(prof@coef[[v]]) + else # then it must be in fixed + newcall$start[[v]] <- as.numeric(prof@call.orig$fixed[v]) + } newobj <- do.call("mle2", newcall) if(class(object) == "fitsad") return (new("fitsad", newobj, sad=object@sad, distr=object@distr, trunc=object@trunc))