Find optimal fits with step-functions having jumps at given candidate positions for all possible subset sizes.
steppath(y, ..., max.blocks)
# S3 method for default
steppath(y, x = 1:length(y), x0 = 2 * x[1] - x[2], max.cand = NULL,
family = c("gauss", "gaussvar", "poisson", "binomial", "gaussKern"), param = NULL,
weights = rep(1, length(y)), cand.radius = 0, ..., max.blocks = max.cand)
# S3 method for stepcand
steppath(y, ..., max.blocks = sum(!is.na(y$number)))
# S3 method for steppath
[[(x, i)
# S3 method for steppath
length(x)
# S3 method for steppath
print(x, ...)
# S3 method for steppath
logLik(object, df = NULL, nobs = object$cand$rightIndex[nrow(object$cand)], ...)For steppath an object of class steppath, i.e. a list with components
pathA list of length length(object) where the ith element contains the best fit by a step-function having i-1 jumps (i.e. i blocks), given by the candidates indices
costA numeric vector of length length(object) giving the value of the cost functional corresponding to the solutions.
candAn object of class stepcand giving the candidates among which the jumps were selected.
[[.steppath returns the fit with i blocks as an object of class stepfit; length.steppath the maximum number of blocks for which a fit has been computed. logLik.stepfit returns an object of class logLik giving the likelihood of the data given the fits corresponding to cost, e.g. for use with AIC.
for steppath:
either an object of class stepcand for steppath.stepcand or a numeric vector containing the serial data for steppath.default
for steppath.default which calls stepcand; see there
single integer giving the maximal number of blocks to find; defaults to number of candidates (note: there will be one block more than the number of jumps
for generic methods only
for methods on a steppath object x or object:
the object
if this is an integer returns the fit with i blocks as an object of class stepcand, else the standard behaviour of a list
the number of estimated parameters: by default the number of blocks for families poisson and binomial, one more (for the variance) for family gauss
the number of observations used for estimating
Friedrich, F., Kempe, A., Liebscher, V., Winkler, G. (2008) Complexity penalized M-estimation: fast computation. Journal of Computational and Graphical Statistics 17(1), 201--224.
# simulate 5 blocks (4 jumps) within a total of 100 data points
b <- c(sort(sample(1:99, 4)), 100)
f <- rep(rnorm(5, 0, 4), c(b[1], diff(b)))
# add Gaussian noise
x <- f + rnorm(100)
# find 10 candidate jumps
cand <- stepcand(x, max.cand = 10)
cand
# compute solution path
path <- steppath(cand)
path
plot(x)
lines(path[[5]], col = "red")
# compare result having 5 blocks with truth
fit <- path[[5]]
fit
logLik(fit)
AIC(logLik(fit))
cbind(fit, trueRightEnd = b, trueLevel = unique(f))
# for poisson observations
y <- rpois(100, exp(f / 10) * 20)
# compute solution path, compare result having 5 blocks with truth
cbind(steppath(y, max.cand = 10, family = "poisson")[[5]],
trueRightEnd = b, trueIntensity = exp(unique(f) / 10) * 20)
# for binomial observations
size <- 10
z <- rbinom(100, size, pnorm(f / 10))
# compute solution path, compare result having 5 blocks with truth
cbind(steppath(z, max.cand = 10, family = "binomial", param = size)[[5]],
trueRightEnd = b, trueIntensity = pnorm(unique(f) / 10))
# an example where stepcand is not optimal but indices found are close to optimal ones
blocks <- c(rep(0, 9), 1, 3, rep(1, 9))
blocks
stepcand(blocks, max.cand = 3)[,c("rightEnd", "value", "number")]
# erroneously puts the "1" into the right block in the first step
steppath(blocks)[[3]][,c("rightEnd", "value")]
# putting the "1" in the middle block is optimal
steppath(blocks, max.cand = 3, cand.radius = 1)[[3]][,c("rightEnd", "value")]
# also looking in the 1-neighbourhood remedies the problem
Run the code above in your browser using DataLab