Nothing
#' Two stage multiple change point detection in linear models.
#'
#' This function provides a two‐stage procedure for simultaneously detecting multiple
#' change‐points in linear models. In the cutting stage, the change‐point problem is
#' converted into a model selection problem so that a modern model selection method
#' can be applied. In the refining stage, the change‐points obtained in the cutting
#' stage are finalized via a refining method. The tuning parameter lambda is chosen
#' by BIC.
#'
#'
#' @param Y an response vector.
#' @param X the n-by-p design matrix.
#' @param method the method to be used by lasso, adaptive lasso, mcp or scad.
#' See \code{\link[plus]{plus}} in R packages \pkg{plus} for details.
#' @param c ceiling(c*sqrt(length(Y))) is the length of each segments in spliting stage.
#' @return tsmcplm returns an object of class "tsmcplm".
#' An object of class "tsmcplm" is a list containing the following components:
#' @return \item{change.points}{estimators of change points.}
#' @export tsmcplm
#' @import plus lars
#' @seealso plus lars
#' @references {Jin, B., Wu, Y., & Shi, X. (2016). Consistent two‐stage multiple
#' change‐point detection in linear models. \emph{Canadian Journal of Statistics},
#' 44(2), 161-179.}
#' @examples ## example 1 : mean shift model
#' ## true change point location:
#' ## 100, 130, 150, 230, 250, 400, 440, 650, 760, 780, 810
#' Y <- rnorm(1000, 0, 0.5) +
#' c(rep(0,100), rep(4,30),rep(-1,20), rep(2,80), rep(-2,20),
#' rep(3,150), rep(-1, 40), rep(1,210), rep(5,110), rep(2,20),
#' rep(7,30), rep(3,190))
#' ts.plot(Y)
#'
#' ##estimate change points
#' tsmcplm(Y = Y, X = NULL, method = "adapt", c = 0.3)
#'
#'
#' ## example 2: linear model:
#' ## a periodic auto correlation series with period 122 and
#' ## order of auto correlation 1
#'
#' ###
#' ## true change point location:
#' ## 200, 350, 450, 550, 700, and 850
#'
#'
#' n=1000
#' y <- rnorm(n)
#' for (t in 2:n) {
#' y[t] <- cos(t*pi/61) + 3*sin(t*pi/61) + 0.5*y[t-1] + (2*sin(t*pi/61) +
#' 0.1 * y[t-1])*(200 < t)+ (2* cos(t*pi/61) - 4 *sin(t*pi/61) -
#' 0.6* y[t-1] )*(350 < t) + (2* sin(t*pi/61) +
#' 0.7* y[t-1] )*(450 < t) + (-3* sin(t*pi/61) -
#' 0.3* y[t-1] )*(550 < t) + (-3* cos(t*pi/61) +
#' 5* sin(t*pi/61))* (700 < t) + (3* cos(t*pi/61) -
#' 5* sin(t*pi/61) - 0.4*y[t-1] )* (850 < t) + rnorm(1)
#'
#' }
#' ts.plot(y)
#'
#' x <- sapply(2:n, function(t){cbind(cos(t*pi/61), sin(t*pi/61), y[t-1])}, simplify = FALSE)
#' x <- do.call(rbind, x)
#' tsmcplm(Y = y[-1], X = x, method = "adapt", c = 2)
#'
#'
tsmcplm <- function(Y, X, method = c("lasso", "adapt", "mcp", "scad"), c) {
# library("lars")
n <- length(Y)
m <- ceiling(c * sqrt(n))
q <- floor(n/m)
K_temp <- matrix(0, nrow = q, ncol = q, byrow = TRUE)
p <- dim(X)[2]
if (length(p) == 0) {
p <- 0
}
########### transform###################
X_temp <- cbind(rep(1, n), X)
Y_temp <- c(Y)
for (i in 1:q) K_temp[i, 1:i] <- rep(1, i)
x <- NULL
y <- NULL
x[[1]] <- as.matrix(X_temp[1:((n - (q - 1) * m)), ])
y[[1]] <- Y_temp[1:((n - (q - 1) * m))]
for (i in 2:q) {
x[[i]] <- as.matrix(X_temp[(n - (q - i + 1) * m + 1):((n - (q -
i) * m)), ])
y[[i]] <- Y_temp[(n - (q - i + 1) * m + 1):((n - (q - i) * m))]
}
X_temp1 <- lapply(1:length(x), function(j, mat, list) kronecker(K_temp[j,
, drop = FALSE], x[[j]]), mat = K_temp, list = x)
Xn <- do.call("rbind", X_temp1)
w1 <- stats::lm(Y_temp ~ Xn - 1)$coefficients
w <- w1
for (i in 1:q) {
w[c(((i - 1) * (p + 1) + 1):(i * (p + 1)))] <- sum(w1[c(((i - 1) *
(p + 1) + 1):(i * (p + 1)))]^2)
}
#################### adlasso#########################
if (method == "adapt") {
X_temp <- scale(Xn, center = FALSE, scale = 1/w) # x1 times the weights with group
object <- lars(X_temp, Y_temp, type = "lasso", intercept = FALSE,
normalize = FALSE)
# get min BIC
bic <- log(n) * object$df + n * log(as.vector(object$RSS)/n) # rss/n version
step.bic2 <- which.min(bic) # step with min BIC
coeff <- coef.lars(object, s = step.bic2, mode = "step")
adcp.coef.s <- sum(abs(coeff))
adcp.coef.v.m <- abs(matrix(c(coeff), q, (p + 1), byrow = T))
adcp.coef.m <- c(apply(adcp.coef.v.m, 1, max))
adcp.cp <- which(adcp.coef.m != 0)
if (length(adcp.cp) > 1) {
for (i in 2:length(adcp.cp)) {
if (adcp.cp[i] - adcp.cp[i - 1] == 1)
adcp.cp[i] <- 0
}
}
adcp.cp1 <- adcp.cp[adcp.cp > 1 & adcp.cp < q]
d1 <- length(adcp.cp1)
if (d1 == 0) {
adcpcss.cp <- 0
}
if (d1 >= 1) {
# step 4 find change points
adcpcss.cp <- NULL
adcp.cp1 <- c(0, adcp.cp1, q + 1)
for (i in 1:d1) {
y1 <- NULL
x1 <- NULL
for (k in (adcp.cp1[i + 1] - 1):(adcp.cp1[i + 1] + 1)) # for(k in (adcp.cp1[i]+1):(adcp.cp1[i+2]-1))
{
y1 <- c(y1, y[[k]])
x1 <- rbind(x1, as.matrix(x[[k]]))
}
# using cusum find change point.
cp <- css(y1, x1)
if (cp == 0)
next
if (cp != 0) {
if (adcp.cp1[i + 1] == 0)
adcpcss.cp <- c(adcpcss.cp, cp)
if (adcp.cp1[i + 1] > 0)
adcpcss.cp <- c(adcpcss.cp, cp + n - (q - adcp.cp1[i +
1] + 2) * m)
# mcpcss.cp=c(mcpcss.cp,cp+n-(q-mcp.cp1[i])*m)
}
}
if (length(adcpcss.cp) == 0)
adcpcss.cp <- 0 else adcpcss.cp <- adcpcss.cp
}
tt <- which(abs(diff(adcpcss.cp)) < 10)
if (length(tt) > 0)
adcpcss.cp <- adcpcss.cp[-tt]
if (length(adcpcss.cp) >= 1 & min(adcpcss.cp) > 0) {
cp1 <- c(0, adcpcss.cp, n)
q <- length(cp1) - 1
K_temp <- matrix(0, nrow = q, ncol = q, byrow = TRUE)
for (i in 1:q) K_temp[i, 1:i] <- rep(1, i)
X_temp <- cbind(rep(1, n), X)
Y_temp <- c(Y)
x <- NULL
y <- NULL
for (t in 2:length(cp1)) {
x[[t - 1]] <- as.matrix(X_temp[(cp1[t - 1] + 1):cp1[t],
])
y[[t - 1]] <- Y_temp[(cp1[t - 1] + 1):cp1[t]]
}
cp0 <- NULL
for (t in 2:(length(cp1) - 1)) {
x1 <- rbind(x[[t - 1]], x[[t]])
y1 <- c(y[[t - 1]], y[[t]])
tt <- css(y1, x1)
if (tt != 0)
cp0 <- c(cp0, tt + cp1[t - 1])
}
adcpcss.cp <- cp0
tt <- which(abs(diff(adcpcss.cp)) < 10)
if (length(tt) > 0)
adcpcss.cp <- adcpcss.cp[-tt]
}
}
################################# lasso#######################
if (method == "lasso") {
X_temp <- Xn # x1 times the weights with group
object <- lars(X_temp, Y_temp, type = "lasso", intercept = FALSE,
normalize = FALSE)
# get min BIC
bic <- log(n) * object$df + n * log(as.vector(object$RSS)/n) # rss/n version
step.bic2 <- which.min(bic) # step with min BIC
coeff <- coef.lars(object, s = step.bic2, mode = "step")
adcp.coef.s <- sum(abs(coeff))
adcp.coef.v.m <- abs(matrix(c(coeff), q, (p + 1), byrow = T))
adcp.coef.m <- c(apply(adcp.coef.v.m, 1, max))
adcp.cp <- which(adcp.coef.m != 0)
if (length(adcp.cp) > 1) {
for (i in 2:length(adcp.cp)) {
if (adcp.cp[i] - adcp.cp[i - 1] == 1)
adcp.cp[i] <- 0
}
}
adcp.cp1 <- adcp.cp[adcp.cp > 1 & adcp.cp < q]
d1 <- length(adcp.cp1)
if (d1 == 0) {
adcpcss.cp <- 0
}
if (d1 >= 1) {
# step 4 find change points
adcpcss.cp <- NULL
adcp.cp1 <- c(0, adcp.cp1, q + 1)
for (i in 1:d1) {
y1 <- NULL
x1 <- NULL
for (k in (adcp.cp1[i + 1] - 1):(adcp.cp1[i + 1] + 1)) # for(k in (adcp.cp1[i]+1):(adcp.cp1[i+2]-1))
{
y1 <- c(y1, y[[k]])
x1 <- rbind(x1, as.matrix(x[[k]]))
}
# using cusum find change point.
cp <- css(y1, x1)
if (cp == 0)
next
if (cp != 0) {
if (adcp.cp1[i + 1] == 0)
adcpcss.cp <- c(adcpcss.cp, cp)
if (adcp.cp1[i + 1] > 0)
adcpcss.cp <- c(adcpcss.cp, cp + n - (q - adcp.cp1[i +
1] + 2) * m)
# mcpcss.cp=c(mcpcss.cp,cp+n-(q-mcp.cp1[i])*m)
}
}
if (length(adcpcss.cp) == 0)
adcpcss.cp <- 0 else adcpcss.cp <- adcpcss.cp
}
tt <- which(abs(diff(adcpcss.cp)) < 10)
if (length(tt) > 0)
adcpcss.cp <- adcpcss.cp[-tt]
if (length(adcpcss.cp) >= 1 & min(adcpcss.cp) > 0) {
cp1 <- c(0, adcpcss.cp, n)
q <- length(cp1) - 1
K_temp <- matrix(0, nrow = q, ncol = q, byrow = TRUE)
for (i in 1:q) K_temp[i, 1:i] <- rep(1, i)
X_temp <- cbind(rep(1, n), X)
Y_temp <- c(Y)
x <- NULL
y <- NULL
for (t in 2:length(cp1)) {
x[[t - 1]] <- as.matrix(X_temp[(cp1[t - 1] + 1):cp1[t],
])
y[[t - 1]] <- Y_temp[(cp1[t - 1] + 1):cp1[t]]
}
cp0 <- NULL
for (t in 2:(length(cp1) - 1)) {
x1 <- rbind(x[[t - 1]], x[[t]])
y1 <- c(y[[t - 1]], y[[t]])
tt <- css(y1, x1)
if (tt != 0)
cp0 <- c(cp0, tt + cp1[t - 1])
}
adcpcss.cp <- cp0
tt <- which(abs(diff(adcpcss.cp)) < 10)
if (length(tt) > 0)
adcpcss.cp <- adcpcss.cp[-tt]
}
}
############################### mcp
if (method == "mcp") {
object <- plus(Xn, Y_temp, method = "mc+", gamma = 2.4, intercept = F,
normalize = F, eps = 1e-30)
# step 1 estimate coef using BIC.
bic <- log(dim(Xn)[1]) * object$dim + dim(Xn)[1] * log(as.vector((1 -
object$r.square) * sum(Y_temp^2))/length(Y_temp))
step.bic <- which.min(bic)
mcp.coef <- coef.plus(object, lam = object$lam[step.bic])
mcp.coef.s <- sum(abs(mcp.coef))
mcp.coef.v.m <- abs(matrix(c(mcp.coef), q, (p + 1), byrow = T))
mcp.coef.m <- c(apply(mcp.coef.v.m, 1, max))
mcp.cp <- which(mcp.coef.m != 0)
if (length(mcp.cp) > 1) {
for (i in 2:length(mcp.cp)) {
if (mcp.cp[i] - mcp.cp[i - 1] == 1)
mcp.cp[i] <- 0
}
}
mcp.cp1 <- mcp.cp[mcp.cp > 1 & mcp.cp < q]
d1 <- length(mcp.cp1)
if (d1 == 0) {
mcpcss.cp <- 0
adcpcss.cp <- mcpcss.cp
}
if (d1 >= 1) {
# step 4 find change points
mcpcss.cp <- NULL
mcp.cp1 <- c(0, mcp.cp1, q + 1)
for (i in 1:d1) {
y1 <- NULL
x1 <- NULL
for (k in (mcp.cp1[i + 1] - 1):(mcp.cp1[i + 1] + 1)) # for(k in (mcp.cp1[i]+1):(mcp.cp1[i+2]-1))
{
y1 <- c(y1, y[[k]])
x1 <- rbind(x1, as.matrix(x[[k]]))
}
# using cusum find change point.
cp <- css(y1, x1)
if (cp == 0)
next
if (cp != 0) {
if (mcp.cp1[i + 1] == 0)
mcpcss.cp <- c(mcpcss.cp, cp)
if (mcp.cp1[i + 1] > 0)
mcpcss.cp <- c(mcpcss.cp, cp + n - (q - mcp.cp1[i +
1] + 2) * m)
# mcpcss.cp=c(mcpcss.cp,cp+n-(q-mcp.cp1[i])*m)
}
}
if (length(mcpcss.cp) == 0)
adcpcss.cp <- 0 else adcpcss.cp <- mcpcss.cp
}
tt <- which(abs(diff(adcpcss.cp)) < 10)
if (length(tt) > 0)
adcpcss.cp <- adcpcss.cp[-tt]
if (length(adcpcss.cp) >= 1 & min(adcpcss.cp) > 0) {
n <- length(Y)
cp1 <- c(0, adcpcss.cp, n)
q <- length(cp1) - 1
K_temp <- matrix(0, nrow = q, ncol = q, byrow = TRUE)
for (i in 1:q) K_temp[i, 1:i] <- rep(1, i)
X_temp <- cbind(rep(1, n), X)
Y_temp <- c(Y)
x <- NULL
y <- NULL
for (t in 2:length(cp1)) {
x[[t - 1]] <- as.matrix(X_temp[(cp1[t - 1] + 1):cp1[t],
])
y[[t - 1]] <- Y_temp[(cp1[t - 1] + 1):cp1[t]]
}
cp0 <- NULL
for (t in 2:(length(cp1) - 1)) {
x1 <- rbind(x[[t - 1]], x[[t]])
y1 <- c(y[[t - 1]], y[[t]])
tt <- css(y1, x1)
if (tt != 0)
cp0 <- c(cp0, tt + cp1[t - 1])
}
adcpcss.cp <- cp0
tt <- which(abs(diff(adcpcss.cp)) < 10)
if (length(tt) > 0)
adcpcss.cp <- adcpcss.cp[-tt]
}
}
############################### scad
if (method == "scad") {
object <- plus(Xn, Y_temp, method = "scad", gamma = 2.4, intercept = F,
normalize = F, eps = 1e-30)
# step 1 estimate coef using BIC.
bic <- log(dim(Xn)[1]) * object$dim + dim(Xn)[1] * log(as.vector((1 -
object$r.square) * sum(Y_temp^2))/length(Y_temp))
step.bic <- which.min(bic)
mcp.coef <- coef.plus(object, lam = object$lam[step.bic])
mcp.coef.s <- sum(abs(mcp.coef))
mcp.coef.v.m <- abs(matrix(c(mcp.coef), q, (p + 1), byrow = T))
mcp.coef.m <- c(apply(mcp.coef.v.m, 1, max))
mcp.cp <- which(mcp.coef.m != 0)
if (length(mcp.cp) > 1) {
for (i in 2:length(mcp.cp)) {
if (mcp.cp[i] - mcp.cp[i - 1] == 1)
mcp.cp[i] <- 0
}
}
mcp.cp1 <- mcp.cp[mcp.cp > 1 & mcp.cp < q]
d1 <- length(mcp.cp1)
if (d1 == 0) {
mcpcss.cp <- 0
adcpcss.cp <- mcpcss.cp
}
if (d1 >= 1) {
# step 4 find change points
mcpcss.cp <- NULL
mcp.cp1 <- c(0, mcp.cp1, q + 1)
for (i in 1:d1) {
y1 <- NULL
x1 <- NULL
for (k in (mcp.cp1[i + 1] - 1):(mcp.cp1[i + 1] + 1)) # for(k in (mcp.cp1[i]+1):(mcp.cp1[i+2]-1))
{
y1 <- c(y1, y[[k]])
x1 <- rbind(x1, as.matrix(x[[k]]))
}
# using cusum find change point.
cp <- css(y1, x1)
if (cp == 0)
next
if (cp != 0) {
if (mcp.cp1[i + 1] == 0)
mcpcss.cp <- c(mcpcss.cp, cp)
if (mcp.cp1[i + 1] > 0)
mcpcss.cp <- c(mcpcss.cp, cp + n - (q - mcp.cp1[i +
1] + 2) * m)
# mcpcss.cp=c(mcpcss.cp,cp+n-(q-mcp.cp1[i])*m)
}
}
if (length(mcpcss.cp) == 0)
adcpcss.cp <- 0 else adcpcss.cp <- mcpcss.cp
}
tt <- which(abs(diff(adcpcss.cp)) < 10)
if (length(tt) > 0)
adcpcss.cp <- adcpcss.cp[-tt]
if (length(adcpcss.cp) >= 1 & min(adcpcss.cp) > 0) {
n <- length(Y)
cp1 <- c(0, adcpcss.cp, n)
q <- length(cp1) - 1
K_temp <- matrix(0, nrow = q, ncol = q, byrow = TRUE)
for (i in 1:q) K_temp[i, 1:i] <- rep(1, i)
X_temp <- cbind(rep(1, n), X)
Y_temp <- c(Y)
x <- NULL
y <- NULL
for (t in 2:length(cp1)) {
x[[t - 1]] <- as.matrix(X_temp[(cp1[t - 1] + 1):cp1[t],
])
y[[t - 1]] <- Y_temp[(cp1[t - 1] + 1):cp1[t]]
}
cp0 <- NULL
for (t in 2:(length(cp1) - 1)) {
x1 <- rbind(x[[t - 1]], x[[t]])
y1 <- c(y[[t - 1]], y[[t]])
tt <- css(y1, x1)
if (tt != 0)
cp0 <- c(cp0, tt + cp1[t - 1])
}
adcpcss.cp <- cp0
tt <- which(abs(diff(adcpcss.cp)) < 10)
if (length(tt) > 0)
adcpcss.cp <- adcpcss.cp[-tt]
}
}
if (length(adcpcss.cp) == 0)
return(0) else return(change.points = adcpcss.cp)
}
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.