R/anotherBkp.R

In mpierrejean/jointseg: Joint Segmentation of Multivariate (Copy Number) Signals

#' Get best candidate change point
#' 
#' Get best candidate change point according to binary segmentation
#' 
#' Contrary to \code{oneBkp}, \code{anotherBkp} handles missing values (NA:s).
#' 
#' @param Y A \code{n*p} matrix, \code{p} signals of length \code{n} to be
#' segmented (centered by column)
#' @param weightFUN A \code{function} returning a \code{(n-1)*1} vector of
#' weights for the candidate change point positions. Default weights yield the
#' likelihood ratio test (LRT) statistic for the identification of a single
#' change point.
#' 
#' @param verbose A \code{logical} value: should extra information be output ?
#' Defaults to \code{FALSE}.
#' @return A numerical value, a candidate breakpoint position
#' @author Morgane Pierre-Jean and Pierre Neuvial
#' @keywords internal
#' @examples
#' 
#' p <- 2
#' n <- 100
#' 
#' sim <- randomProfile(n, 1, 1, p)
#' Y <- sim$profile
#' bkp <- anotherBkp(Y)
#' print(bkp)
#' print(oneBkp(Y))
#' ##  stopifnot(identical(oneBkp(Y), bkp))
#' plotSeg(Y, list(sim$bkp, bkp))
#' 
#' ## robustness to NA:s
#' h <- 2
#' idxs <- seq(from=max(sim$bkp[1]-h, 1), min(sim$bkp[1]+h, n))
#' Y[idxs, p] <- NA
#' oneBkp(Y)  ## does not work
#' bkp <- anotherBkp(Y)  ## works
#' bkp-sim$bkp
#' 
#' @export anotherBkp
anotherBkp <- function(Y, weightFUN=defaultWeights, verbose=FALSE) {
    if (!is.matrix(Y)) {
        stop("Y is not a matrix, please check dimension of Y")
    }
    C <- apply(Y, 2, getUnivStat, weightFUN=weightFUN)
    cNorm <- rowSums(C^2, na.rm=TRUE)
    which.max(cNorm)
}

############################################################################
## HISTORY:
## 2013-01-30
## o Now an internal function (not exported anymore).
## o Added 'jointSeg:::' to example because function is not exported anymore.
## 2013-01-09
## o Replaced 'jump' by 'bkp'.
## 2012-12-30
## o Created from 'oneBkp'.
############################################################################