Nothing
stepFit <- function(y, q = NULL, alpha = NULL, x = 1:length(y), x0 = 2 * x[1] - x[2], family = NULL, intervalSystem = NULL,
lengths = NULL, confband = FALSE, jumpint = confband, ...) {
if (missing(y)) {
stop("argument 'y' must be given")
}
.RemoveAdditionalArgsPF <- function(family, y, n, ..., penalty, alpha, stat, r, weights, options,
seed, rand.gen, messages)
.parametricFamily(family = family, y = y, n = n, ...)
data <- .RemoveAdditionalArgsPF(family = family, y = y, n = length(y), ...)
intervalSystem <- .intervalSystem(intervalSystem = intervalSystem, lengths = lengths, data = data)
if (length(x) != data$n) {
stop("x and y must have the same length")
}
if (!is.numeric(x) || !all(is.finite(x))) {
stop("x must be a finite numeric vector")
}
if (any(x[-1] - x[-data$n] <= 0)) {
stop("x must be an increasing sequence")
}
if (!is.numeric(x0) || length(x0) != 1 || !is.finite(x0) || x0 >= x[1]) {
stop("x0 must be a single finite numeric smaller than x[1]")
}
.RemoveAdditionalArgsCV <- function(q, alpha, data, output, intervalSystem, ..., nq,
sd, covariances, correlations, filter)
.critVal(q = q, alpha = alpha, data = data, output = output, intervalSystem = intervalSystem, ...)
q <- .RemoveAdditionalArgsCV(q = q, alpha = alpha, data = data, output = "vector",
intervalSystem = intervalSystem, ...)
criticalValues <- rep(Inf, data$n)
criticalValues[intervalSystem$lengths] <- q
if (length(confband) != 1 || !is.logical(confband) || is.na(confband)) {
stop("confband must be a single logical (not NA)")
}
if (length(jumpint) != 1 || !is.logical(jumpint) || is.na(jumpint)) {
stop("jumpint must be a single logical (not NA)")
}
if (confband) {
jumpint <- TRUE
ret <- .callRoutines(observations = data$y, routineType = 5L,
argumentsListRoutine = list(q = criticalValues),
dataType = data$type, argumentsListData = data$argumentsList,
intervalSystemType = intervalSystem$type,
argumentsListIntervalSystem = intervalSystem$argumentsList)
} else {
if (jumpint) {
ret <- .callRoutines(observations = data$y, routineType = 4L,
argumentsListRoutine = list(q = criticalValues),
dataType = data$type, argumentsListData = data$argumentsList,
intervalSystemType = intervalSystem$type,
argumentsListIntervalSystem = intervalSystem$argumentsList)
} else {
if (intervalSystem$type < 10L) {
ret <- .callRoutines(observations = data$y, routineType = 3L,
argumentsListRoutine = list(q = criticalValues),
dataType = data$type, argumentsListData = data$argumentsList,
intervalSystemType = intervalSystem$type,
argumentsListIntervalSystem = intervalSystem$argumentsList)
} else {
ret <- .callRoutines(observations = data$y, routineType = 4L,
argumentsListRoutine = list(q = criticalValues),
dataType = data$type, argumentsListData = data$argumentsList,
intervalSystemType = intervalSystem$type,
argumentsListIntervalSystem = intervalSystem$argumentsList)
}
}
}
fit <- stepfit(cost = attr(ret, "cost"), family = data$family, value = ret$value,
param = data$argumentsListData,
x0 = x0, leftEnd = x[ret$leftIndex], rightEnd = x[ret$rightIndex],
leftIndex = ret$leftIndex, rightIndex = ret$rightIndex)
if (jumpint) {
fit$leftEndLeftBound <- x[c(1L, ret$leftConfInt)]
fit$leftEndRightBound <- x[c(1L, ret$rightConfInt)]
fit$rightEndLeftBound <- x[c(ret$leftConfInt - 1L, data$n)]
fit$rightEndRightBound <- x[c(ret$rightConfInt - 1L, data$n)]
fit$leftIndexLeftBound <- c(1L, ret$leftConfInt)
fit$leftIndexRightBound <- c(1L, ret$rightConfInt)
fit$rightIndexLeftBound <- c(ret$leftConfInt - 1L, data$n)
fit$rightIndexRightBound <- c(ret$rightConfInt - 1L, data$n)
}
if (confband) {
band <- data.frame(x = x, lower = ret$lowerBand, upper = ret$upperBand)
attr(band, "x0") <- x0
class(band) <- c("confband", class(band))
attr(fit, "confband") <- band
}
fit
}
"stepfit" <-
function(cost, family, value, param = NULL, leftEnd, rightEnd, x0, leftIndex = leftEnd, rightIndex = rightEnd)
{
ret <- stepblock(value, leftEnd, rightEnd, x0)
ret$leftIndex <- leftIndex
ret$rightIndex <- rightIndex
attr(ret, "cost") <- cost
attr(ret, "family") <- family
attr(ret, "param") <- param
class(ret) <- c("stepfit", class(ret))
ret
}
"print.stepfit" <-
function(x, ...)
{
cat("\n")
cat("Fitted step function of family", attr(x, "family"), "containing", nrow(x), "blocks\n\n")
cat("domain: (", attr(x, "x0"), ",", x$rightEnd[nrow(x)], "]\n")
cat("range: [", min(x$value), ",", max(x$value), "]\n")
cat("cost:", attr(x, "cost"), "\n")
if(!is.null(attr(x, "param"))) {
cat("parameter:\n")
print(attr(x, "param"))
}
cat("\n")
}
"[.stepfit" <-
function (x, i, j, drop = if(missing(i)) TRUE else if(missing(j)) FALSE else length(j) == 1, refit = FALSE)
{
ret <- NextMethod("[.")
# refit
if(!identical(refit, FALSE)) {
if(missing(i)) i <- 1:nrow(x)
ret$value <- switch(attr(x, "family"),
gauss = diff(c(0, ret$cumSum)) / diff(c(attr(ret, "x0"), ret$rightIndex)),
gaussKern = if(nrow(ret) == 1 ) diff(c(0, ret$cumSum)) / diff(c(attr(ret, "x0"), ret$rightIndex)) else {
param <- attr(x, "param")
r <- ret$rightIndex
ir <- r[-length(r)] # inner positions
n <- r[length(r)] # last position
kl <- length(param$kern)
kj <- param$jump
s <- param$step
d <- c(ir - kj, n) - c(0, ir + kl - kj) + c(rep(sum((1 - s)^2), length(r) - 1), 0) + c(0, rep(sum(s^2), length(r) - 1))
ld <- length(d)
sd <- sum( s * ( 1 - s ) )
XX <- diag(d)
XX[cbind(1:(ld-1), 2:ld)] <- sd
XX[cbind(2:ld, 1:(ld-1))] <- sd
dif <- diff(c(0, ret$rightIndex))
close <- min(dif) <= kl
Xy <- c(0, ret$lXy[-nrow(ret)]) + ret$rcXy - c(0, ret$lcXy[-nrow(ret)]) + ret$rXy
if(identical(refit, TRUE) | !close) {
if(close) warning("jumps closer than filter length")
} else {
ss <- c(0, s, 1)
revss <- c(1, rev(s), 0)
for(i in which(dif < kl)) {
if(i > 1 & i < n) {
# compute design matrix for neighbouring blocks
neigh <- which(r + kj > r[i-1] & c(0, ir) <= r[i] + kl - kj )
neighn <- length(neigh)
neighi <- max(r[neigh[1]] - kj + 1, 1):min(r[neigh[neighn] - 1] + kl - kj, n)
tX <- outer(1:neighn, neighi, function(k,ind) {
# print(data.frame(k = k, neighk = neigh[k], rneighk = r[neigh[k]], ind = ind, test = ind <= r[neigh[k]] - kj, s = pmin(pmax(c(0,r)[neigh[k]] + kl - kj + 1 - ind, 0), kl + 1) + 1, revss = revss[pmin(pmax(c(0,r)[neigh[k]] + kl - kj + 1 - ind, 0), kl + 1) + 1], s1 = pmin(pmax(ind + kj - r[neigh[k]], 0), kl + 1) + 1, ss = ss[pmin(pmax(ind + kj - r[neigh[k]], 0), kl + 1) + 1]))
revss[pmin(pmax(c(0,r)[neigh[k]] + kl - kj + 1 - ind, 0), kl + 1) + 1] - ss[pmin(pmax(ind + kj - r[neigh[k]], 0), kl + 1) + 1]
})
# print(tX)
iXX <- tX %*% t(tX)
# apply to off-diagonal elements including or crossing i
XX[neigh[neigh <= i],neigh[neigh >= i]] <- iXX[which(neigh <= i),which(neigh >= i)]
XX[neigh[neigh >= i],neigh[neigh <= i]] <- iXX[which(neigh >= i),which(neigh <= i)]
Xy[i] <- tX[which(neigh == i),] %*% refit[neighi]
}
}
}
ret$value <- as.numeric(solve(XX, Xy))
},
poisson = diff(c(0, ret$cumSum)) / diff(c(attr(ret, "x0"), ret$rightIndex)),
binomial = diff(c(0, ret$cumSum)) / diff(c(attr(ret, "x0"), ret$rightIndex)) / attr(ret, "param")
)
class(ret) <- class(x)
}
ret
}
"plot.stepfit" <-
function(x, dataspace = TRUE, ...)
{
if(attr(x, "family") == "binomial" && dataspace) {
x$value <- x$value * attr(x, "param")
}
NextMethod()
}
"lines.stepfit" <-
function(x, dataspace = TRUE, ...)
{
if(attr(x, "family") == "binomial" && dataspace) {
x$value <- x$value * attr(x, "param")
}
NextMethod()
}
"fitted.stepfit" <-
function(object, ...)
{
if(attr(object, "family") == "gaussKern" && nrow(object) > 1) {
k <- attr(object, "param")$kern
l <- length(k)
j <- attr(object, "param")$jump
ret <- rep(object$value, object$rightIndex - object$leftIndex + 1 + c(l - j - 1, rep(0, length(object$rightIndex) - 2), j))
convolve(ret, rev(k), conj = TRUE, type = "filter")
} else if(attr(object, "family") == "binomial") {
rep(object$value * attr(object, "param"), object$rightIndex - object$leftIndex + 1)
} else {
rep(object$value, object$rightIndex - object$leftIndex + 1)
}
}
"residuals.stepfit" <-
function(object, y, ...)
{
fit <- fitted(object)
if(length(fit) != length(y)) stop("data and fit differ in length")
switch(attr(object, "family"),
binomial = {
y - attr(object, "param") * fit
},
gaussvar = {
ifelse(fit == 0, ifelse(y^2 == 0, 0, Inf), log(y^2 / fit)) # sort of residuals: log( y^2 / sigma^2 )
},
{
y - fit
}
)
}
"logLik.stepfit" <-
function(object, df = NULL, nobs = object$rightIndex[nrow(object)], ...)
{
family <- switch(attr(object, "family"),
gaussKern = "gauss",
# gaussInhibit = "gauss",
# gaussInhibitBoth = "gauss",
attr(object, "family")
)
ret <- switch(family,
gauss = {
if(is.null(df)) df <- nrow(object) + 1 # variance has also been estimated
-nobs / 2 * (1 + log(2 * pi * attr(object, "cost") / nobs))
},
{
if(is.null(df)) df <- nrow(object)
attr(object, "cost")
}
)
attr(ret, "df") <- df
attr(ret, "nobs") <- nobs
class(ret) <- "logLik"
ret
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.