steppath {stepR}

Solution path of step-functions

Description

Find optimal fits with step-functions having jumps at given candidate positions for all possible subset sizes.

Usage

steppath(y, ..., max.blocks)
## Default S3 method:
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 class 'stepcand'
steppath(y, ..., max.blocks = sum(!is.na(y$number)))
## S3 method for class 'steppath'
x[[i]]
## S3 method for class 'steppath'
length(x)
## S3 method for class 'steppath'
print(x, ...)
## S3 method for class 'steppath'
logLik(object, df = NULL, nobs = object$cand$rightIndex[nrow(object$cand)], ...)

Arguments

for steppath:

y

either an object of class stepcand for steppath.stepcand or a numeric vector containing the serial data for steppath.default

x, x0, max.cand, family, param, weights, cand.radius

for steppath.default which calls stepcand; see there

max.blocks

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:

object

the object

i

if this is an integer returns the fit with i blocks as an object of class stepcand, else the standard behaviour of a list

df

the number of estimated parameters: by default the number of blocks for families poisson and binomial, one more (for the variance) for family gauss

nobs

the number of observations used for estimating

Value

For steppath an object of class steppath, i.e. a list with components

path

A 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

cost

A numeric vector of length length(object) giving the value of the cost functional corresponding to the solutions.

cand

An 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.

References

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.

See Also

stepcand, stepfit, family, logLik, AIC

Examples

# 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