# R/profile_mlePP.r In NHPoisson: Modelling and Validation of Non Homogeneous Poisson Processes

setMethod("profile",
signature(fitted = "mlePP"),
function (fitted, which = 1L:p, maxsteps = 100,
alpha = 0.01, zmax = sqrt(qchisq(1 - alpha, 1L)), del = zmax/5,
trace = FALSE, ...)
{

onestep <- function(step) {
bi <- B0[i] + sgn * step * del * std.err[i]
fix <- list(bi)
names(fix) <- pi
call$fixed <- c(fix, fix0) pfit <- tryCatch(eval.parent(call, 2L), error = identity) if (inherits(pfit, "error")) return(NA) else { zz <- 2 * (pfit@min - fitted@min) ri <- pv0 ri[, names(pfit@coef)] <- pfit@coef ri[, pi] <- bi if (zz > -0.001) zz <- max(zz, 0) else stop("profiling has found a better solution, so original fit had not converged") z <- sgn * sqrt(zz) pvi <<- rbind(pvi, ri) zi <<- c(zi, z) } if (trace) cat(bi, z, "\n") z } summ <- summary(fitted) std.err <- summ@coef[, "Std. Error"] Pnames <- names(B0 <- fitted@coef) pv0 <- t(as.matrix(B0)) p <- length(Pnames) prof <- vector("list", length = length(which)) names(prof) <- Pnames[which] call <- fitted@call call$modSim <- TRUE
call$modCI <- FALSE call$dplot <- FALSE
ndeps <- eval.parent(call$control$ndeps)
parscale <- eval.parent(call$control$parscale)
fix0 <- eval.parent(call$fixed) 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 lastz <- 0 while ((step <- step + 1) < maxsteps && abs(z) < zmax) { z <- onestep(step) 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) if (is.na(z) || abs(z) > zmax) break } } else if (length(zi) < 5L) { mxstep <- step - 1L step <- 0.5 while ((step <- step + 1L) < mxstep) onestep(step) } } 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)
}
)


## Try the NHPoisson package in your browser

Any scripts or data that you put into this service are public.

NHPoisson documentation built on Feb. 19, 2020, 5:07 p.m.