R/sampleEntropy.R In nonlinearTseries: Nonlinear Time Series Analysis

#' Sample Entropy (also known as Kolgomorov-Sinai Entropy)
#' @description
#' The Sample Entropy measures the complexity of a time series. Large values of
#' the Sample Entropy indicate high complexity whereas that smaller values
#' characterize more regular signals.
#' @details  The sample entropy is computed using:
#' \deqn{h_q(m,r) = log(C_q(m,r)/C_{q}(m+1,r))}{hq(m,r) = log(Cq(m,r)/Cq(m+1,r)),}
#' where \emph{m} is the embedding dimension and \emph{r} is the radius of the
#' neighbourhood. When computing the correlation dimensions we use the linear
#' regions from the correlation sums in order to do the estimates. Similarly,
#' the sample entropy \eqn{h_q(m,r)}{hq(m,r)} should not change for both
#' various \emph{m} and \emph{r}.
#' @param corrDim.object A \emph{corrDim} object from which the Sample Entropy
#' of the time series characterized by \emph{corrDim} shall be estimated.
#' @param do.plot do.plot Logical value. If TRUE (default value), a plot of the
#' sample entropy is shown.
#' @param ... Additional plotting arguments.
#' @return A \emph{sampleEntropy} object that contains a list storing the
#' sample entropy (\emph{sample.entropy}), the embedding dimensions
#' entropy has been computed, and the order of the sample entropy
#' (\emph{entr.order}). The sample entropy is stored as a matrix in which
#' each row contains the computations for a given embedding dimension and
#' each column stores the computations for a given radius.
#' @references H. Kantz  and T. Schreiber: Nonlinear Time series Analysis
#' (Cambridge university press)
#' @examples
#' \dontrun{
#' x=henon(n.sample = 15000, n.transient = 100, a = 1.4, b = 0.3,
#'         start = c(0.78,0.8165), do.plot = FALSE)$x #' #' cd=corrDim(time.series=x, #' min.embedding.dim=2,max.embedding.dim=9, #' corr.order=2,time.lag=1, #' min.radius=0.05,max.radius=1, #' n.points.radius=100, #' theiler.window=20, #' do.plot=TRUE) #' #' use.col = c("#999999", "#E69F00", "#56B4E9", "#009E73", #' "#F0E442", "#0072B2", "#D55E00", "#CC79A7") #' se=sampleEntropy(cd,do.plot=TRUE,col=use.col, #' type="l",xlim=c(0.1,1), #' add.legend=T) #' se.est = estimate(se, #' regression.range = c(0.4,0.6), #' use.embeddings = 6:9,col=use.col,type="b") #' print(se.est) #' cat("Expected K2 = ",0.325," Estimated = ",mean(se.est),"\n") #' } #' @author Constantino A. Garcia #' @rdname sampleEntropy #' @export sampleEntropy #' @exportClass sampleEntropy sampleEntropy = function (corrDim.object, do.plot = TRUE, ...) { if (!inherits(corrDim.object, "corrDim")) { stop("corrDim.object should be of class corrDim") } radius = radius(corrDim.object) corr.matrix = corrMatrix(corrDim.object) embeddings = embeddingDims(corrDim.object) number.embeddings = length(embeddings) - 1 entropy = matrix(0, nrow = number.embeddings, ncol = length(radius)) for (i in 1:number.embeddings) { entropy[i,] = log(corr.matrix[i, ] / corr.matrix[i + 1, ]) } dimnames(entropy) = list(head(embeddings, -1), radius) sample.entropy = list(sample.entropy = entropy, embedding.dims = head(embeddings, -1), entr.order = nlOrder(corrDim.object), radius = radius ) class(sample.entropy) = "sampleEntropy" attr(sample.entropy,"time.lag") = attr(corrDim.object,"time.lag") attr(sample.entropy,"id") = attr(corrDim.object,"id") attr(sample.entropy,"theiler.window") = attr(corrDim.object,"theiler.window") if (do.plot) { plot(sample.entropy, ...) } sample.entropy } #' Returns the sample entropy function \eqn{h_q(m,r)} used for the computations #' of the sample entropy. #' @param x A \emph{sampleEntropy} object. #' @return A numeric matrix representing the sample entropy function #' \eqn{h_q(m,r)} obtained in #' the sample entropy computations represented #' by the \emph{sampleEntropy} object. #' @seealso \code{\link{sampleEntropy}} #' @export sampleEntropyFunction sampleEntropyFunction = function(x){ UseMethod("sampleEntropyFunction") } #' @return The \emph{sampleEntropyFunction} returns the sample entropy function #' \eqn{h_q(m,r)} used for the computations. The sample #' entropy function is represented by a matrix. Each row represents a given #' embedding dimension whereas that each column representes a different radius. #' @rdname sampleEntropy #' @export sampleEntropyFunction.sampleEntropy = function(x){ x$sample.entropy
}

#' @return The \emph{nlOrder} function returns the order of the sample entropy.
#' @rdname sampleEntropy
#' @export
nlOrder.sampleEntropy = function(x){
x$entr.order } #' @return The \emph{radius} function returns the radius on which the sample #' entropy function has been evaluated. #' @rdname sampleEntropy #' @export radius.sampleEntropy = function(x){ radius.default(x) } #' @return The \emph{embeddingDims} function returns the embedding dimensions #' on which the sample entropy function has been evaluated. #' @rdname sampleEntropy #' @export embeddingDims.sampleEntropy = function(x){ embeddingDims.default(x) } #' @return The \emph{plot} function shows the graphics for the sample entropy. #' @rdname sampleEntropy #' @param main A title for the plot. #' @param xlab A title for the x axis. #' @param ylab A title for the y axis. #' @param type Type of plot (see \code{\link[graphics]{plot}}). #' @param col Vector of colors for each of the dimensions of the plot. #' @param pch Vector of symbols for each of the dimensions of the plot #' @param ylim Numeric vector of length 2, giving the y coordinates range.. #' @param log A character string which contains "x" if the x axis is to be #' logarithmic, "y" if the y axis is to be logarithmic and "xy" or "yx" if #' both axes are to be logarithmic. #' @param add.legend add a legend to the plot? #' @export plot.sampleEntropy = function(x, main = NULL, xlab = NULL, ylab = NULL, type = "l", col=NULL, pch=NULL, ylim=NULL, log = "x", add.legend = T, ...) { if (is.null(xlab)) { xlab = expression(epsilon) } if (is.null(ylab)) { ylab = expression(h[q]*"("*epsilon*")") } if (is.null(main)) { main = bquote( "Sample entropy (q = " * .(x$entr.order) * ")    " * h[.(x$entr.order)] * "(" * epsilon * ")" ) } number.embeddings = length(x$embedding.dims)

# obtain vector of graphical parameters if not specified
col = vectorizePar(col,number.embeddings)
pch = vectorizePar(pch,number.embeddings)

on.exit(par(current.par))
layout(rbind(1, 2), heights = c(8, 2))
}
for (i in 1:number.embeddings) {
if (i == 1) {
if (is.null(ylim)) {
ylim = range(x$sample.entropy) } plot(x$radius, x$sample.entropy[1, ], xlab = xlab, ylab = ylab, main = main, type = type, pch = pch[[i]], col = col[[i]], ylim = ylim, log = log, ...) } else { lines(x$radius, x$sample.entropy[i, ], type = type, pch = pch[[i]], col = col[[i]], ...) } } if (add.legend) { par(mar = c(0, 0, 0, 0)) plot.new() legend("center", "groups", ncol = ceiling(number.embeddings / 2), col = col, pch = pch, lty = rep(1, number.embeddings), lwd = rep(2.5, number.embeddings), bty = "n", legend = x$embedding.dims, title = "Embedding dimension")
}

}

#' @details For each embedding dimension the sample
#' entropy is estimated by averaging  \deqn{h_q(m,r) = log(C_q(m,r)/C_{q}(m+1,r))}{hq(m,r) = log(Cq(m,r)/Cq(m+1,r))}
#' over the range specified by \emph{regression range} in the \emph{estimate}
#' function.
#' @param x A \emph{sampleEntropy} object.
#' @param regression.range Vector with 2 components denoting the range where the
#'  function will perform linear regression.
#' @param use.embeddings A numeric vector specifying which embedding dimensions
#' should the \emph{estimate} function use to compute
#' the sample entropy.
#' @param fit.lty The type of line to plot the regression lines.
#' @param fit.lwd	The width of the line for the regression lines.
#' @param fit.col A vector of colors to plot the regression lines.
#' @return The  \emph{estimate} function returns a vector storing the sample
#' entropy estimate for each embedding dimension.
#' @rdname sampleEntropy
#' @export
estimate.sampleEntropy = function(x, regression.range = NULL, do.plot = TRUE,
use.embeddings = NULL, fit.col = NULL,
fit.lty = 2, fit.lwd = 2,
if (is.null(regression.range)) {
regression.range = range(x$radius) } if (is.null(use.embeddings)) { use.embeddings = x$embedding.dims
}
#select only embeddings used in the object
use.embeddings = intersect(as.numeric(use.embeddings), x$embedding.dims) len.use.embeddings = length(use.embeddings) if (len.use.embeddings < 1) { stop("The embeddings specified are not present!\n") } #range indx = which(x$radius >= regression.range[[1]] &
x$radius <= regression.range[[2]]) if (len.use.embeddings == 1) { sample.entropy.estimate = mean( x$sample.entropy[as.character(use.embeddings), indx]
)
} else {
sample.entropy.estimate = apply(
x$sample.entropy[as.character(use.embeddings), indx], MARGIN = 1, FUN = mean) } if (do.plot) { #create new sample entropy object with the embedding.dims and radius range used x = list(sample.entropy = x$sample.entropy[as.character(use.embeddings), ],
embedding.dims = use.embeddings, entr.order = x$entr.order, radius = x$radius)
class(x) = "sampleEntropy"
plotSampleEntropyEstimate(x, sample.entropy.estimate,
fit.col = fit.col, fit.lty = fit.lty, fit.lwd = fit.lwd,
}
sample.entropy.estimate
}

plotSampleEntropyEstimate = function(sampleEntropy.object,
sample.entropy.estimate, main = NULL,
xlab = NULL, ylab = NULL, type = "l",
pch = NULL, ylim = NULL, col = NULL,
fit.col, fit.lty, fit.lwd,
number.embeddings = length(sampleEntropy.object$embedding.dims) # obtain vector of graphical parameters if not specified col = vectorizePar(col,number.embeddings) pch = vectorizePar(pch,number.embeddings) fit.col = vectorizePar(fit.col,number.embeddings,col) if (add.legend) { current.par = par(no.readonly = TRUE) on.exit(par(current.par)) layout(rbind(1, 2), heights = c(8, 2)) } plot(sampleEntropy.object, main = main, xlab = xlab, ylab = ylab, type = type, pch = pch, ylim = ylim, col = col, add.legend = F, ...) for (i in 1:number.embeddings) { if (i == 1) { abline(h = sample.entropy.estimate[[i]], col = fit.col[[i]], lty = fit.lty, lwd = fit.lwd) } else { abline(h = sample.entropy.estimate[[i]], col = fit.col[[i]], lty = fit.lty, lwd = fit.lwd) } } if (add.legend) { par(mar = c(0, 0, 0, 0)) plot.new() legend("center", "groups", ncol = ceiling(number.embeddings / 2), col = col, pch = pch, lty = rep(1, number.embeddings), lwd = rep(2.5, number.embeddings), bty = "n", legend = sampleEntropy.object$embedding.dims,