Nothing
R/wvar.R
In wv: Wavelet Variance
Defines functions
unitConversion
wvar
wvar.lts
wvar.gts
wvar.ts
wvar.imu
wvar.default
create_wvar
print.wvar
summary.wvar
plot.wvar
Documented in
create_wvar
plot.wvar
print.wvar
summary.wvar
unitConversion
wvar
wvar.default
wvar.gts
wvar.imu
wvar.lts
wvar.ts
#' @title Convert Unit of Time Series Data
#' @description Manipulate the units of time to different ones
#' @keywords internal
#' @param x A \code{vector} containing the values on x-axis.
#' @param from.unit A \code{string} indicating the unit which the data is converted from.
#' @param to.unit A \code{string} indicating the unit which the data is converted to.
#' @details
#' The supported units are "ns"(nanosecond), "ms"(millisecond), "sec", "min", "hour", "day", "month", and "year".
#' Make sure \code{from.unit} and \code{to.unit} are not \code{NULL} before it is passed to this function.
#' @return A \code{list} with the following structure:
#' \itemize{
#' \item "x": Data
#' \item "converted": A \code{boolean} indicating whether conversion is made
#' }
#' @export
#' @examples
#' x = seq(60, 3600, 60)
#' unitConversion(x, 'sec', 'min')
#' y = 1:10
#' unitConversion(y, 'hour', 'sec')
unitConversion = function(x, from.unit, to.unit){
#ns, ms, second, min, hour, day, month, year
unit = c(ns = 1, ms = 2,se = 3, mi = 4, ho = 5, da = 6, mo = 7, ye = 8)
#assume 1 month = 30 days
ratio = c(1E6, 1E3, 60, 60, 24, 30, 12)
from.unit.1 = substr(from.unit, 1, 2)
to.unit.1 = substr(to.unit, 1, 2)
#check unit:
no.convert = F
if(from.unit.1 == to.unit.1){no.convert = T}
if(is.na(unit[from.unit.1]) ) {
message = paste('No such unit: ', from.unit, '. Supported units are "ns"(nanosecond), "ms"(millisecond), "sec", "min", "hour", "day", "month", and "year". Conversion is terminated.', sep = '')
warning(message); no.convert = T}
if(is.na(unit[to.unit.1]) ) {
message = paste('No such unit: ', to.unit, '. Supported units are "ns"(nanosecond), "ms"(millisecond), "sec", "min", "hour", "day", "month", and "year". Conversion is terminated.', sep = '')
warning(message); no.convert = T}
if(!no.convert){
#print out warning when day is convert to month, or month is converted to day.
conversionRange = unit[from.unit.1] : unit[to.unit.1]
if(6 %in% conversionRange && 7 %in% conversionRange){
warning('Unit conversion might be wrong because this function simply assumes 1 month = 30 days.')
}
}
if(!no.convert){
if(unit[from.unit.1] > unit[to.unit.1]){
temp = ratio[unit[to.unit.1]: (unit[from.unit.1]-1)]
multiplier = prod(temp)
x = x*multiplier
}else{
temp = ratio[unit[from.unit.1]: (unit[to.unit.1]-1) ]
multiplier = prod(temp)
x = x/multiplier
}
}
obj = list(x = x, converted = !no.convert)
return(obj)
}
#' @title Wavelet Variance
#' @description Calculates the (MO)DWT wavelet variance
#' @param x A \code{vector} with dimensions N x 1.
#' @param decomp A \code{string} that indicates whether to use a "dwt" or "modwt" decomposition.
#' @param filter A \code{string} that specifies which wavelet filter to use.
#' @param nlevels An \code{integer} that indicates the level of decomposition. It must be less than or equal to floor(log2(length(x))).
#' @param alpha A \code{double} that specifies the significance level which in turn specifies the \eqn{1-\alpha} confidence level.
#' @param robust A \code{boolean} that triggers the use of the robust estimate.
#' @param eff A \code{double} that indicates the efficiency as it relates to an MLE.
#' @param freq A \code{numeric} that provides the rate of samples.
#' @param from.unit A \code{string} indicating the unit from which the data is converted.
#' @param to.unit A \code{string} indicating the unit to which the data is converted.
#' @param ... Further arguments passed to or from other methods.
#' @return A \code{list} with the structure:
#' \itemize{
#' \item "variance": Wavelet Variance
#' \item "ci_low": Lower CI
#' \item "ci_high": Upper CI
#' \item "robust": Robust active
#' \item "eff": Efficiency level for Robust calculation
#' \item "alpha": p value used for CI
#' \item "unit": String representation of the unit
#' }
#' @details
#' The default value of \code{nlevels} will be set to \eqn{\left\lfloor {{{\log }_2}\left( {length\left( x \right)} \right)} \right\rfloor}{floor(log2(length(x)))}, unless otherwise specified.
#' @author James Balamuta, Justin Lee and Stephane Guerrier
#' @rdname wvar
#' @examples
#' set.seed(999)
#' x = rnorm(100)
#'
#' # Default
#' wvar(x)
#'
#' # Robust
#' wvar(x, robust = TRUE, eff=0.3)
#'
#' # Classical
#' wvar(x, robust = FALSE, eff=0.3)
#'
#' # 90% Confidence Interval
#' wvar(x, alpha = 0.10)
#' @export
wvar = function(x, ...) {
if (sum(class(x) %in% "gts") == 1){
x = as.numeric(x)
}
UseMethod("wvar")
}
#' @rdname wvar
#' @export
wvar.lts = function(x, decomp = "modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, to.unit = NULL, ...){
warning('`lts` object is detected. This function can only operate on the combined process.')
freq = attr(x, 'freq')
unit = attr(x, 'unit')
x = x[,ncol(x)]
wvar.default(x, decomp, filter, nlevels, alpha, robust, eff, freq = freq, from.unit = unit, to.unit = to.unit)
}
#' @rdname wvar
#' @export
wvar.gts = function(x, decomp="modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, to.unit = NULL, ...){
freq = attr(x, 'freq')
unit = attr(x, 'unit')
x = x[,1]
wvar.default(x, decomp, filter, nlevels, alpha, robust, eff, freq = freq, from.unit = unit, to.unit = to.unit)
}
#' @rdname wvar
#' @export
wvar.ts = function(x, decomp="modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, to.unit = NULL, ...){
freq = attr(x, 'tsp')[3]
unit = NULL
wvar.default(x, decomp, filter, nlevels, alpha, robust, eff, freq = freq, from.unit = unit, to.unit = to.unit)
}
#' @rdname wvar
#' @export
wvar.imu = function(x, decomp="modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, to.unit = NULL, ...){
# Retrive sensor name
if (!is.null(attr(x, "stype"))){
sensor_name = attr(x, "stype")
}else{
warning("Unknown sensor name. IMU object is missing some information.")
sensor_name = NULL
}
# Retrive freq
if (!is.null(attr(x, "freq"))){
freq = attr(x, "freq")
}else{
warning("Unknown frequency. IMU object is missing some information. Freq is set to 1 by default.")
freq = 1
}
# Retrive sample size
if (!is.null(attr(x, "dim"))){
n = attr(x, "dim")[1]
}else{
warning("Unknown sample size. IMU object is missing some information.")
n = NULL
}
# Retrive col names
if (!is.null(attr(x, "dimnames")[[2]])){
col_names = attr(x, "dimnames")[[2]]
}else{
stop("Unknown colunms names. IMU object is missing some information.")
col_names = NULL
}
# Retrive sensor
if (!is.null(attr(x, "sensor"))){
sensor = attr(x, "sensor")
}else{
warning("Unknown sensor. IMU object is missing some information.")
sensor = NULL
}
# Retrive axis
if (!is.null(attr(x, "axis"))){
ax = attr(x, "axis")
}else{
warning("Unknown axes. IMU object is missing some information.")
ax = NULL
}
# Compute wvar
m = length(col_names)
wvariance = list()
for (i in 1:m){
wvariance[[i]] = wvar.default(x[,i], decomp, filter, nlevels, alpha, robust, eff, freq = freq, to.unit = to.unit)
}
names(wvariance) = col_names
out = list(sensor = sensor_name, freq = freq, n = n, type = sensor, axis = ax, wvar = wvariance)
class(out) = "imu_wvar"
invisible(out)
}
#' @rdname wvar
#' @export
#' @importFrom methods is
wvar.default = function(x, decomp = "modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, freq = 1, from.unit = NULL, to.unit = NULL, ...){
if(is.null(x)){
stop("`x` must contain a value")
}else if((is.data.frame(x) || is.matrix(x))){
if(ncol(x) > 1) stop("There must be only one column of data supplied.")
}
if(decomp == "modwt" && is.null(nlevels)){
nlevels = floor(log2(length(x))-1)
}else if(decomp == "dwt" && is.null(nlevels)){
nlevels = floor(log2(length(x)))
}
# Check Freq
if(!is(freq,"numeric") || length(freq) != 1){ stop("'freq' must be one numeric number.") }
if(freq <= 0) { stop("'freq' must be larger than 0.") }
# Check Unit
all.units = c('ns', 'ms', 'sec', 'second', 'min', 'minute', 'hour', 'day', 'mon', 'month', 'year')
if( (!is.null(from.unit) && !from.unit %in% all.units) || (!is.null(to.unit) && !to.unit %in% all.units) ){
stop('The supported units are "ns", "ms", "sec", "min", "hour", "day", "month", "year". ')
}
if(robust) {
if(eff > 0.99) {
stop("The efficiency specified is too close to the classical case. Use `robust = FALSE`")
}
}
obj = modwt_wvar_cpp(signal=x, nlevels=nlevels, robust=robust, eff=eff, alpha=alpha,
ci_type="eta3", strWavelet=filter, decomp = decomp)
# nlevels may be changed during modwt
nlevels = nrow(obj)
scales = scales_cpp(nlevels)/freq
# NO Unit Conversion
if( is.null(from.unit) && is.null(to.unit)==F ){
warning("'from.unit' is NULL. Unit conversion was not done.")
}
# Unit Conversion
if (!is.null(from.unit)){
if (!is.null(to.unit)){
convert.obj = unitConversion(scales, from.unit = from.unit, to.unit = to.unit)
if (convert.obj$converted) {
# YES Unit Conversion
scales = convert.obj$x
message(paste0('Unit of object is converted from ', from.unit, ' to ', to.unit), appendLF = T)
}
}
}
if(!is.null(from.unit) && !is.null(to.unit)){
unit = to.unit
}else{
unit = from.unit}
create_wvar(obj, decomp, filter, robust, eff, alpha, scales, unit)
}
#' @title Create a \code{wvar} object
#' @description Structures elements into a \code{wvar} object
#' @param obj A \code{matrix} with dimensions N x 3 that contains Wavelet Variance, Lower CI, and Upper CI.
#' @param decomp A \code{string} that indicates whether to use a "dwt" or "modwt" decomposition.
#' @param filter A \code{string} that specifies the type of wavelet filter used in the decomposition.
#' @param robust A \code{boolean} that triggers the use of the robust estimate.
#' @param eff A \code{double} that indicates the efficiency as it relates to an MLE.
#' @param alpha A \code{double} that specifies the significance level which in turn specifies the \eqn{1-\alpha} confidence level.
#' @param scales A \code{vec} that contains the amount of decomposition performed at each level.
#' @param unit A \code{string} that indicates the unit expression of the frequency.
#' @return A \code{list} with the structure:
#' \itemize{
#' \item "variance": Wavelet variance
#' \item "ci_low": Lower CI
#' \item "ci_high": Upper CI
#' \item "robust": Robust active
#' \item "eff": Efficiency level for robust calculation
#' \item "alpha": p value used for CI
#' \item "unit": String representation of the unit
#' }
#' @keywords internal
create_wvar = function(obj, decomp, filter, robust, eff, alpha, scales, unit){
structure(list(variance = obj[,1],
ci_low = obj[,2],
ci_high = obj[,3],
robust = robust,
eff = eff,
alpha = alpha,
scales = scales,
decomp = decomp,
unit = unit,
filter = filter), class = "wvar")
}
#' @title Print Wavelet Variances
#' @description Displays the summary table of wavelet variance.
#' @author James Balamuta
#' @method print wvar
#' @export
#' @keywords internal
#' @param x A \code{wvar} object.
#' @param ... Further arguments passed to or from other methods.
#' @return Summary table
#' @examples
#' set.seed(999)
#' x = rnorm(100)
#' out = wvar(x)
#' print( out )
print.wvar = function(x, ...){
mat = matrix(unlist(x[1:3]),ncol=3,byrow=F)
colnames(mat) = c("Variance", "Low CI", "High CI")
rownames(mat) = x$scales
print(mat)
}
#' @title Summary of Wavelet Variances
#' @description Displays the summary table of wavelet variance accounting for CI values and supplied efficiency.
#' @method summary wvar
#' @export
#' @keywords internal
#' @param object A \code{wvar} object.
#' @param ... Additional arguments affecting the summary produced.
#' @return Summary table and other properties of the object.
#' @author James Balamuta
#' @examples
#' set.seed(999)
#' x = rnorm(100)
#' ret = wvar(x)
#' summary(ret)
summary.wvar = function(object, ...){
name = if(object$robust){
"robust"
}else{
"classical"
}
cat("Results of the wavelet variance calculation using the ",name, " method.\n",sep="")
if(object$robust){
cat("Robust was created using efficiency=",object$eff,"\n",sep="")
}
cat("The confidence interval was generated using (1-",object$alpha,")*100 \n",sep="")
print(object)
}
#' @title Plot Wavelet Variance
#' @description Displays a plot of wavelet variance accounting for CI values and supplied efficiency.
#' @method plot wvar
#' @keywords internal
#' @param x A \code{wvar} object.
#' @param units A \code{string} that specifies the units of time plotted on the x axis.
#' @param xlab A \code{string} that gives a title for the x axis.
#' @param ylab A \code{string} that gives a title for the y axis.
#' @param main A \code{string} that gives an overall title for the plot.
#' @param col_wv A \code{string} that specifies the color of the wavelet variance line.
#' @param col_ci A \code{string} that specifies the color of the confidence interval polygon.
#' @param nb_ticks_x An \code{integer} that specifies the maximum number of ticks for the x-axis.
#' @param nb_ticks_y An \code{integer} that specifies the maximum number of ticks for the y-axis.
#' @param legend_position A \code{string} that specifies the position of the legend (use \code{legend_position = NA} to remove legend).
#' @param ci_wv A \code{boolean} that determines whether a confidence interval polygon will be drawn.
#' @param point_cex A \code{double} that specifies the size of each symbol to be plotted.
#' @param point_pch A \code{double} that specifies the symbol type to be plotted.
#' @param ... Additional arguments affecting the plot.
#' @return Plot of wavelet variance and confidence interval for each scale.
#' @author Stephane Guerrier, Nathanael Claussen, and Justin Lee
#' @export
#' @examples
#' set.seed(999)
#' n = 10^4
#' Xt = rnorm(n)
#' wv = wvar(Xt)
#' plot(wv)
#' plot(wv, main = "Simulated white noise", xlab = "Scales")
#' plot(wv, units = "sec", legend_position = "topright")
#' plot(wv, col_wv = "darkred", col_ci = "pink")
plot.wvar = function(x, units = NULL, xlab = NULL, ylab = NULL, main = NULL,
col_wv = NULL, col_ci = NULL, nb_ticks_x = NULL, nb_ticks_y = NULL,
legend_position = NULL, ci_wv = NULL, point_cex = NULL,
point_pch = NULL, ...){
# Labels
if (is.null(xlab)){
if (is.null(units)){
xlab = expression(paste("Scale ", tau, sep =""))
}else{
xlab = bquote(paste("Scale ", tau, " [", .(units), "]", sep = " "))
}
}
if (is.null(ylab)){
ylab = expression(paste("Wavelet Variance ", nu^2, sep = ""))
}else{
ylab = ylab
}
# Main Title
if (is.null(main)){
main = "Haar Wavelet Variance Representation"
}
# Line and CI colors
if (is.null(col_wv)){
col_wv = "darkblue"
}
if (is.null(col_ci)){
col_ci = hcl(h = 210, l = 65, c = 100, alpha = 0.2)
}
# Range
x_range = range(x$scales)
x_low = floor(log2(x_range[1]))
x_high = ceiling(log2(x_range[2]))
y_range = range(c(x$ci_low, x$ci_high))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
# Axes
if (is.null(nb_ticks_x)){
nb_ticks_x = 6
}
if (is.null(nb_ticks_y)){
nb_ticks_y = 5
}
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > nb_ticks_x){
x_ticks = x_low + ceiling((x_high - x_low)/(nb_ticks_x + 1))*(0:nb_ticks_x)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(2^ .(i))))
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels <- sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
# Legend Position
if (is.null(legend_position)){
if (which.min(abs(c(y_low, y_high) - log2(x$variance[1]))) == 1){
legend_position = "topleft"
}else{
legend_position = "bottomleft"
}
}
# Main Plot
plot(NA, xlim = x_range, ylim = y_range, xlab = xlab, ylab = ylab,
log = "xy", xaxt = 'n', yaxt = 'n', bty = "n", ann = FALSE)
win_dim = par("usr")
par(new = TRUE)
plot(NA, xlim = x_range, ylim = 10^c(win_dim[3], win_dim[4] + 0.09*(win_dim[4] - win_dim[3])),
xlab = xlab, ylab = ylab, log = "xy", xaxt = 'n', yaxt = 'n', bty = "n")
win_dim = par("usr")
# Add Grid
abline(v = 2^x_ticks, lty = 1, col = "grey95")
abline(h = 10^y_ticks, lty = 1, col = "grey95")
# Add Title
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
polygon(x_vec, y_vec, col = "grey95", border = NA)
text(x = 10^mean(c(win_dim[1], win_dim[2])), y = 10^(win_dim[4] - 0.09/2*(win_dim[4] - win_dim[3])), main)
# Add Axes and Box
lines(x_vec[1:2], rep(10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])),2), col = 1)
#y_ticks = y_ticks[(2^y_ticks) < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))]
y_labels = y_labels[1:length(y_ticks)]
box()
axis(1, at = 2^x_ticks, labels = x_labels, padj = 0.3)
axis(2, at = 10^y_ticks, labels = y_labels, padj = -0.2)
# CI for WV
if (ci_wv == TRUE || is.null(ci_wv)){
polygon(c(x$scales, rev(x$scales)), c(x$ci_low, rev(x$ci_high)),
border = NA, col = col_ci)
}
# Add legend
CI_conf = 1 - x$alpha
if (x$robust == TRUE){
wv_title_part1 = "Empirical Robust WV "
}else{
wv_title_part1 = "Empirical WV "
}
if (!is.na(legend_position)){
if (legend_position == "topleft"){
legend_position = 10^c(1.1*win_dim[1], 0.98*(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])))
legend(x = legend_position[1], y = legend_position[2],
legend = c(as.expression(bquote(paste(.(wv_title_part1), hat(nu)^2))),
as.expression(bquote(paste("CI(",hat(nu)^2,", ",.(CI_conf),")")))),
pch = c(16, 15), lty = c(1, NA), col = c(col_wv, col_ci), cex = 1, pt.cex = c(1.25, 3), bty = "n")
}else{
if (legend_position == "topright"){
legend_position = 10^c(0.7*win_dim[2], 0.98*(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])))
legend(x = legend_position[1], y = legend_position[2],
legend = c(as.expression(bquote(paste(.(wv_title_part1), hat(nu)^2))),
as.expression(bquote(paste("CI(",hat(nu)^2,", ",.(CI_conf),")")))),
pch = c(16, 15), lty = c(1, NA), col = c(col_wv, col_ci), cex = 1, pt.cex = c(1.25, 3), bty = "n")
}else{
legend(legend_position,
legend = c(as.expression(bquote(paste(.(wv_title_part1), hat(nu)^2))),
as.expression(bquote(paste("CI(",hat(nu)^2,", ",.(CI_conf),")")))),
pch = c(16, 15), lty = c(1, NA), col = c(col_wv, col_ci), cex = 1, pt.cex = c(1.25, 3), bty = "n")
}
}
}
# Add WV
lines(x$scales, x$variance, type = "l", col = col_wv, pch = 16)
if (is.null(point_pch)){
point_pch = 16
}
if (is.null(point_cex)){
point_cex = 1.25
}
lines(x$scales, x$variance, type = "p", col = col_wv, pch = point_pch, cex = point_cex)
}
#' @title Plot Wavelet Variance based on IMU Data
#' @description Displays a plot of wavelet variance accounting for CI values and supplied efficiency.
#' @method plot imu_wvar
#' @keywords internal
#' @param x A \code{wvar} object.
#' @param xlab A \code{string} that gives a title for the x axis.
#' @param ylab A \code{string} that gives a title for the y axis.
#' @param main A \code{string} that gives an overall title for the plot.
#' @param col_wv A \code{string} that specifies the color of the wavelet variance line.
#' @param col_ci A \code{string} that specifies the color of the confidence interval polygon.
#' @param nb_ticks_x An \code{integer} that specifies the maximum number of ticks for the x-axis.
#' @param nb_ticks_y An \code{integer} that specifies the maximum number of ticks for the y-axis.
#' @param ci_wv A \code{boolean} that determines whether a confidence interval polygon will be drawn.
#' @param point_cex A \code{double} that specifies the size of each symbol to be plotted.
#' @param point_pch A \code{double} that specifies the symbol type to be plotted.
#' @param ... Additional arguments affecting the plot.
#' @return Plot of wavelet variance and confidence interval for each scale.
#' @author Stephane Guerrier and Yuming Zhang
#' @export
#' @examples
#' data("kvh1750_wv")
#' plot(kvh1750_wv)
plot.imu_wvar = function(x, xlab = NULL, ylab = NULL, main = NULL,
col_wv = NULL, col_ci = NULL, nb_ticks_x = NULL, nb_ticks_y = NULL,
ci_wv = NULL, point_cex = NULL, point_pch = NULL, ...){
type = unique(x$type)
if ("Gyroscope" %in% type){
gyro_index = which(x$type == "Gyroscope")
}else{
gyro_index = NULL
}
if ("Accelerometer" %in% type){
accel_index = which(x$type == "Accelerometer")
}else{
accel_index = NULL
}
ncol = length(unique(x$axis))
nrow = length(type)
m = length(x$wvar)
J = length(x$wvar[[1]]$variance)
# remove negative CI values
index_to_remove = c()
for (i in 1:m) {
if(length(which(x$wvar[[i]]$lci<0)) > 0){
index_to_remove = c(index_to_remove, which(x$wvar[[i]]$lci<0))
}
}
if (!is.null(index_to_remove)){
index_to_remove = unique(index_to_remove)
index_to_keep = which(seq(1:J) != index_to_remove)
}else{
index_to_keep = 1:J
}
J = length(index_to_keep)
scales = x$wvar[[1]]$scales[index_to_keep]
ci_up = ci_lw = av = matrix(NA, J, m)
for (i in 1:m){
ci_up[,i] = x$wvar[[i]]$ci_high[index_to_keep]
ci_lw[,i] = x$wvar[[i]]$ci_low[index_to_keep]
av[,i] = x$wvar[[i]]$variance[index_to_keep]
}
# Axes
if (is.null(nb_ticks_x)){
nb_ticks_x = 6
}
if (is.null(nb_ticks_y)){
nb_ticks_y = 5
}
# Range
x_range = range(scales)
x_low = floor(log10(x_range[1]))
x_high = ceiling(log10(x_range[2]))
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > nb_ticks_x){
x_ticks = x_low + ceiling((x_high - x_low)/(nb_ticks_x + 1))*(0:nb_ticks_x)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(10^ .(i))))
# Line and CI colors
if (is.null(col_wv)){
col_wv = "darkblue"
}
if (is.null(col_ci)){
col_ci = hcl(h = 210, l = 65, c = 100, alpha = 0.2)
}
if (is.null(point_pch)){
point_pch = 16
}
if (is.null(point_cex)){
point_cex = 1.25
}
# Main Title
if (is.null(main)){
main = paste("Wavelet Variance Representation - ", x$sensor, " @ ", x$freq, " Hz", sep="")
}
# Labels
if (is.null(xlab)){
xlab = bquote(paste("Averaging time ", tau, " [sec]", sep = " "))
}
if (is.null(ylab)){
ylab = expression(paste("Wavelet Variance ", nu, sep = ""))
}
# Main plot
par(omi=rep(1, 4), mar=c(0,0,0,0), mfrow=c(nrow,ncol))
# Gyro
if (!is.null(gyro_index)){
y_range = c(min(ci_lw[,gyro_index]), max(ci_up[,gyro_index]))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels <- sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
for (i in seq_along(gyro_index)){
plot(NA, xlim = range(scales), ylim = y_range, xaxt="n", yaxt="n", log = "xy", bty = "n")
box(col = "grey")
mtext(paste("Axis - ", x$axis[gyro_index][i], sep = ""), 3, line = 0.5)
if (i == 1){
axis(2, at = 10^y_ticks, labels = y_labels, padj = -0.2, cex = 1.25)
}
if (i == 1){
mtext("Gyroscope", 2, line = 4.5)
mtext(ylab, 2, line = 2.5)
}
abline(h = 10^y_ticks, col = "grey85")
abline(v = 10^x_ticks, col = "grey85")
# CI for AD
if(ci_wv == TRUE || is.null(ci_wv)){
polygon(c(scales, rev(scales)), c(ci_lw[,gyro_index[i]], rev(ci_up[,gyro_index[i]])),
border = NA, col = col_ci)
}
# Add AD
lines(scales, (av[,gyro_index[i]]), type = "l", col = col_wv, pch = 16)
lines(scales, (av[,gyro_index[i]]), type = "p", col = col_wv, pch = point_pch, cex = point_cex)
if (is.null(accel_index)){
axis(1, at = 10^x_ticks, labels = x_labels, padj = -0.2, cex = 1.25)
}
}
}
# Accel
if (!is.null(accel_index)){
y_range = c(min(ci_lw[,accel_index]), max(ci_up[,accel_index]))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels <- sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
for (i in seq_along(accel_index)){
plot(NA, xlim = range(scales), ylim = y_range, xaxt="n", yaxt="n", log = "xy", bty = "n")
box(col = "grey")
if (i == 1){
axis(2, at = 10^y_ticks, labels = y_labels, padj = -0.2, cex = 1.25)
}
if (i == 1){
mtext("Accelerometer", 2, line = 4.5)
mtext(ylab, 2, line = 2.5)
}
if (length(accel_index) == 3 && i == 2){
mtext(xlab, 1, line = 3.5)
}
if (is.null(gyro_index)){
mtext(paste("Axis - ", x$axis[gyro_index][i], sep = ""), 3, line = 0.5)
}
abline(h = 10^y_ticks, col = "grey85")
abline(v = 10^x_ticks, col = "grey85")
# CI for AD
if(ci_wv == TRUE || is.null(ci_wv)){
polygon(c(scales, rev(scales)), c(ci_lw[,accel_index[i]], rev(ci_up[,accel_index[i]])),
border = NA, col = col_ci)
}
# Add AD
lines(scales, (av[,accel_index[i]]), type = "l", col = col_wv, pch = 16)
lines(scales, (av[,accel_index[i]]), type = "p", col = col_wv, pch = point_pch, cex = point_cex)
axis(1, at = 10^x_ticks, labels = x_labels, padj = -0.2, cex = 1.25)
}
}
# Add main title
mtext(main, side = 3, line = 3, outer = TRUE)
par(mfrow = c(1,1))
}
#' @title Comparison between classical and robust Wavelet Variances
#' @description Displays a plot of the wavelet variances (classical and robust) for a given time series accounting for CI values.
#' @param x A time series objects.
#' @param eff An \code{integer} that specifies the efficiency of the robust estimator.
#' @param units A \code{string} that specifies the units of time plotted on the x axis.
#' @param xlab A \code{string} that gives a title for the x axis.
#' @param ylab A \code{string} that gives a title for the y axis.
#' @param main A \code{string} that gives an overall title for the plot.
#' @param col_wv A \code{string} that specifies the color of the wavelet variance line.
#' @param col_ci A \code{string} that specifies the color of the confidence interval shade.
#' @param nb_ticks_x An \code{integer} that specifies the maximum number of ticks for the x-axis.
#' @param nb_ticks_y An \code{integer} that specifies the maximum number of ticks for the y-axis.
#' @param legend_position A \code{string} that specifies the position of the legend (use \code{legend_position = NA} to remove legend).
#' @param ... Additional arguments affecting the plot.
#' @return Plot of wavelet variance and confidence interval for each scale.
#' @author Stephane Guerrier, Nathanael Claussen, and Justin Lee
#' @examples
#' set.seed(999)
#' n = 10^4
#' Xt = rnorm(n)
#' wv = wvar(Xt)
#'
#' plot(wv)
#' plot(wv, main = "Simulated white noise", xlab = "Scales")
#' plot(wv, units = "sec", legend_position = "topright")
#' plot(wv, col_wv = "darkred", col_ci = "pink")
#' @export
robust_eda = function(x, eff = 0.6, units = NULL, xlab = NULL, ylab = NULL, main = NULL,
col_wv = NULL, col_ci = NULL, nb_ticks_x = NULL, nb_ticks_y = NULL,
legend_position = NULL, ...){
wv_cl = wvar(x)
wv_rob = wvar(x, robust = TRUE, eff = eff)
# Labels
if (is.null(xlab)){
if (is.null(units)){
xlab = expression(paste("Scale ", tau, sep =""))
}else{
xlab = bquote(paste("Scale ", "", tau , " [", .(units), "]", sep = ""))
}
}
if (is.null(ylab)){
ylab = expression(paste("Wavelet Variance ", nu^2, sep = ""))
}else{
ylab = ylab
}
# Main Title
if (is.null(main)){
main = "Classical vs Robust WV"
}
# Line and CI colors
if (is.null(col_wv)){
col_wv = c("darkblue", "darkorange2")
}
if (is.null(col_ci)){
col_ci = c(hcl(h = 210, l = 65, c = 100, alpha = 0.2), hcl(h = 60, l = 65, c = 100, alpha = 0.2))
}
# Range
x_range = range(wv_cl$scales)
x_low = floor(log2(x_range[1]))
x_high = ceiling(log2(x_range[2]))
y_range = range(c(wv_cl$ci_low, wv_cl$ci_high, wv_rob$ci_low, wv_rob$ci_high))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
# Axes
if (is.null(nb_ticks_x)){
nb_ticks_x = 6
}
if (is.null(nb_ticks_y)){
nb_ticks_y = 5
}
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > nb_ticks_x){
x_ticks = x_low + ceiling((x_high - x_low)/(nb_ticks_x + 1))*(0:nb_ticks_x)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(2^ .(i))))
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels <- sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
# Legend position
if (is.null(legend_position)){
if (which.min(abs(c(y_low, y_high) - log2(wv_rob$variance[1]))) == 1){
legend_position = "topleft"
}else{
legend_position = "bottomleft"
}
}
# Main plot
plot(NA, xlim = x_range, ylim = y_range, xlab = xlab, ylab = ylab,
log = "xy", xaxt = 'n', yaxt = 'n', bty = "n", ann = FALSE)
win_dim = par("usr")
par(new = TRUE)
plot(NA, xlim = x_range, ylim = 10^c(win_dim[3], win_dim[4] + 0.09*(win_dim[4] - win_dim[3])),
xlab = xlab, ylab = ylab, log = "xy", xaxt = 'n', yaxt = 'n', bty = "n")
win_dim = par("usr")
# Add grid
abline(v = 2^x_ticks, lty = 1, col = "grey95")
abline(h = 10^y_ticks, lty = 1, col = "grey95")
# Add title
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
polygon(x_vec, y_vec, col = "grey95", border = NA)
text(x = 10^mean(c(win_dim[1], win_dim[2])), y = 10^(win_dim[4] - 0.09/2*(win_dim[4] - win_dim[3])), main)
# Add Axes and Box
lines(x_vec[1:2], rep(10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])),2), col = 1)
y_ticks = y_ticks[(2^y_ticks) < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))]
y_labels = y_labels[1:length(y_ticks)]
box()
axis(1, at = 2^x_ticks, labels = x_labels, padj = 0.3)
axis(2, at = 10^y_ticks, labels = y_labels, padj = -0.2)
# CI for WV
polygon(c(wv_cl$scales, rev(wv_cl$scales)), c(wv_cl$ci_low, rev(wv_cl$ci_high)),
border = NA, col = col_ci[1])
polygon(c(wv_rob$scales, rev(wv_rob$scales)), c(wv_rob$ci_low, rev(wv_rob$ci_high)),
border = NA, col = col_ci[2])
# Legend Position
if (!is.na(legend_position)){
if (legend_position == "topleft"){
legend_position = 10^c(1.1*win_dim[1], 0.98*(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])))
legend(x = legend_position[1], y = legend_position[2],
legend = c("Classical WV", "Classical CI", "Robust WV", "Robust CI"),
pch = c(16, 15, 16, 15), lty = c(1, NA, 1, NA), col = c(col_wv[1], col_ci[1], col_wv[2], col_ci[2]),
cex = 1, pt.cex = c(1.25, 3, 1.25, 3), bty = "n")
}else{
if (legend_position == "topright"){
legend_position = 10^c(0.7*win_dim[2], 0.98*(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])))
legend(x = legend_position[1], y = legend_position[2],
legend = c("Classical WV", "Classical CI", "Robust WV", "Robust CI"),
pch = c(16, 15, 16, 15), lty = c(1, NA, 1, NA), col = c(col_wv[1], col_ci[1], col_wv[2], col_ci[2]),
cex = 1, pt.cex = c(1.25, 3, 1.25, 3), bty = "n")
}else{
legend(legend_position,
legend = c("Classical WV", "Classical CI", "Robust WV", "Robust CI"),
pch = c(16, 15, 16, 15), lty = c(1, NA, 1, NA), col = c(col_wv[1], col_ci[1], col_wv[2], col_ci[2]),
cex = 1, pt.cex = c(1.25, 3, 1.25, 3), bty = "n")
}
}
}
lines(wv_cl$scales, wv_cl$variance, type = "l", col = col_wv[1], pch = 16)
lines(wv_cl$scales, wv_cl$variance, type = "p", col = col_wv[1], pch = 16, cex = 1.25)
lines(wv_cl$scales, wv_rob$variance, type = "l", col = col_wv[2], pch = 16)
lines(wv_cl$scales, wv_rob$variance, type = "p", col = col_wv[2], pch = 16, cex = 1.25)
}
#' @title Multi-Plot Comparison Between Multiple Wavelet Variances
#' @description
#' This is a helper function for the \code{compare_var()} function.
#' This method accepts the same set of arguments as \code{compare_wvar} and returns a comparision
#' of multiple wavelet variances of different time series accounting for CI values as a set of different plots.
#'
#' @param graph_details List of inputs
#'
#' @author Stephane Guerrier, Justin Lee, and Nathanael Claussen
#' @export
#'
compare_wvar_split = function(graph_details){
old_pars = par(mfrow = c(graph_details$obj_len, graph_details$obj_len),
mar = c(0.5,0.5,0.5,1.5), oma = c(4,4,4,4))
on.exit(par(old_pars))
for (i in 1:graph_details$obj_len){
for (j in 1:graph_details$obj_len){
# Main plot
plot(NA, xlim = graph_details$x_range, ylim = graph_details$y_range, log = "xy", xaxt = 'n',
yaxt = 'n', bty = "n", ann = FALSE)
win_dim = par("usr")
kill_y_tick = graph_details$y_at < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
# Add grid
abline(v = graph_details$x_at, lty = 1, col = "grey95")
abline(h = graph_details$y_at, lty = 1, col = "grey95")
# Add axes and box
box(col = "grey")
# Corner left piece
if (j == 1){
axis(2, at = graph_details$y_at[kill_y_tick],
labels = graph_details$y_labels[kill_y_tick], padj = -0.2, cex.axis = 1/log(graph_details$obj_len))
}
# Corner bottom
if (i == graph_details$obj_len){
axis(1, at = graph_details$x_at, labels = graph_details$x_labels,
padj = 0.1, cex.axis = 1/log(graph_details$obj_len)) # figure out how to size these things for smaller plots
}
# Diag graph
if (i == j){
scales = graph_details$obj_list[[i]]$scales
ci_low = graph_details$obj_list[[i]]$ci_low
ci_high = graph_details$obj_list[[i]]$ci_high
variance = graph_details$obj_list[[i]]$variance
if(graph_details$ci_wv[i] == TRUE){
polygon(c(scales, rev(scales)), c(ci_low, rev(ci_high)),
border = NA, col = graph_details$col_ci[i])
}
lines(scales, variance, type = "l", col = graph_details$col_wv[i], pch = 16)
lines(scales, variance, type = "p", col = graph_details$col_wv[i],
pch = 17, cex = graph_details$point_cex[i]/1.25)
win_dim = par("usr")
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
box()
# if (graph_details$add_legend){
# if (i == j){
# legend(graph_details$legend_position, graph_details$names, bty = "n",
# lwd = 1, pt.cex = graph_details$point_cex, pch = graph_details$point_pch,
# col = graph_details$col_wv, cex=0.7)
# }
# }
}
if (i != j){
scales = graph_details$obj_list[[i]]$scales
ci_low = graph_details$obj_list[[i]]$ci_low
ci_high = graph_details$obj_list[[i]]$ci_high
variance = graph_details$obj_list[[i]]$variance
win_dim = par("usr")
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
if (is.null(graph_details$main[i,j])){
main = paste("WV:", graph_details$names[i],
"vs", graph_details$names[j])
}
box()
if (i < j && graph_details$ci_wv[i] == TRUE){
polygon(c(scales, rev(scales)), c(ci_low, rev(ci_high)),
border = NA, col = graph_details$col_ci[i])
}
lines(scales, variance, type = "l", col = graph_details$col_wv[i], pch = 16)
lines(scales, variance, type = "p", col = graph_details$col_wv[i],
pch = 17, cex = graph_details$point_cex[i]/1.25)
scales = graph_details$obj_list[[j]]$scales
ci_low = graph_details$obj_list[[j]]$ci_low
ci_high = graph_details$obj_list[[j]]$ci_high
variance = graph_details$obj_list[[j]]$variance
if (i < j && graph_details$ci_wv[i] == TRUE){ # don't show confidence intervals
polygon(c(scales, rev(scales)), c(ci_low, rev(ci_high)),
border = NA, col = graph_details$col_ci[j])
}
lines(scales, variance, type = "l", col = graph_details$col_wv[j], pch = 16)
lines(scales, variance, type = "p", col = graph_details$col_wv[j],
pch = 17, cex = graph_details$point_cex[j]/1.25)
}
# Add Details
# @todo: expand win_dim and position $names
if(j==4){
x_vec = 13^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
#mtext(graph_details$names[i], side = 4, line = 0.1, cex = 0.8)
par(xpd = TRUE) #Draw outside plot area
text(x = x_vec[2], y = 0.03, graph_details$names[i], srt = 270, cex = 1.3, col = graph_details$col_wv[i])
par(xpd = FALSE)
}
if(i==1){
mtext(graph_details$names[j], side = 3, line = 0.2, cex = 0.8, col = graph_details$col_wv[j])
}
}
}
mtext(graph_details$ylab, side = 2, line = 2.80, cex = graph_details$cex_labels, outer = T)
mtext(graph_details$xlab, side = 1, line = 2.75, cex = graph_details$cex_labels, outer = T)
if (is.null(graph_details$main)){
main = "Haar Wavelet Variance Representation"
}else{
mtext(main, side = 3, line = 1)
}
}
#' @title Combined Plot Comparison Between Multiple Wavelet Variances
#' @description
#' This is a helper function for the \code{compare_var()} function.
#' This method accepts the same set of arguments as \code{compare_wvar} and returns a single plot
#' that compares multiple wavelet variances of different time series accounting for CI values.
#'
#' @param graph_details List of inputs
#'
#' @author Stephane Guerrier, Justin Lee, and Nathanael Claussen
#' @export
#'
compare_wvar_no_split = function(graph_details){
# Main plot
plot(NA, xlim = graph_details$x_range, ylim = graph_details$y_range, log = "xy", xaxt = 'n',
yaxt = 'n', bty = "n", ann = FALSE)
win_dim = par("usr")
# Main Plot
par(new = TRUE)
plot(NA, xlim = graph_details$x_range, ylim = 10^c(win_dim[3], win_dim[4] + 0.09*(win_dim[4] - win_dim[3])),
log = "xy", xaxt = 'n', yaxt = 'n', bty = "n",
xlab = graph_details$xlab, ylab = graph_details$ylab,
cex.lab = graph_details$cex_labels)
win_dim = par("usr")
# Add Grid
abline(v = graph_details$x_at, lty = 1, col = "grey95")
abline(h = graph_details$y_at, lty = 1, col = "grey95")
# Add Title
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
polygon(x_vec, y_vec, col = "grey95", border = NA)
text(x = 10^mean(c(win_dim[1], win_dim[2])), y = 10^(win_dim[4] - 0.09/2*(win_dim[4] - win_dim[3])),
graph_details$main)
# Add Axes and Box
lines(x_vec[1:2], rep(10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])),2), col = 1)
y_ticks = graph_details$y_ticks[(10^graph_details$y_ticks) < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))]
kill_y_tick = graph_details$y_at < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
box()
axis(1, at = graph_details$x_at, labels = graph_details$x_labels, padj = 0.3)
axis(2, at = graph_details$y_at[kill_y_tick],
labels = graph_details$y_labels[kill_y_tick], padj = -0.2)
for (i in 1:graph_details$obj_len){
scales = graph_details$obj_list[[i]]$scales
ci_low = graph_details$obj_list[[i]]$ci_low
ci_high = graph_details$obj_list[[i]]$ci_high
variance = graph_details$obj_list[[i]]$variance
if (graph_details$ci_wv[i]){
polygon(c(scales, rev(scales)), c(ci_low, rev(ci_high)),
border = NA, col = graph_details$col_ci[i])
}
lines(scales, variance, type = "l", col = graph_details$col_wv[i], pch = 16)
lines(scales, variance, type = "p", col = graph_details$col_wv[i],
pch = graph_details$point_pch[i], cex = graph_details$point_cex[i])
}
if (graph_details$add_legend){
legend(graph_details$legend_position, graph_details$names, bty = "n",
lwd = 1, pt.cex = graph_details$point_cex, pch = graph_details$point_pch,
col = graph_details$col_wv)
}
}
#' @title Comparison Between Multiple Wavelet Variances
#' @description
#' Displays plots of multiple wavelet variances of different time series accounting for CI values.
#'
#' @param ... One or more time series objects.
#' @param split A \code{boolean} that, if TRUE, arranges the plots into a matrix-like format.
#' @param add_legend A \code{boolean} that, if TRUE, adds a legend to the plot.
#' @param units A \code{string} that specifies the units of time plotted on the x axes. Note: This argument will not be used if xlab is specified.
#' @param xlab A \code{string} that gives a title for the x axes.
#' @param ylab A \code{string} that gives a title for the y axes.
#' @param main A \code{string} that gives an overall title for the plot.
#' @param col_wv A \code{string} that specifies the color of the wavelet variance lines.
#' @param col_ci A \code{string} that specifies the color of the confidence interval shade.
#' @param nb_ticks_x An \code{integer} that specifies the maximum number of ticks for the x-axis.
#' @param nb_ticks_y An \code{integer} that specifies the maximum number of ticks for the y-axis.
#' @param legend_position A \code{string} that specifies the position of the legend (use \code{legend_position = NA} to remove legend).
#' @param ci_wv A \code{boolean} that determines whether confidence interval polygons will be drawn.
#' @param point_cex A \code{double} that specifies the size of each symbol to be plotted.
#' @param point_pch A \code{double} that specifies the symbol type to be plotted.
#' @param names A \code{string} that specifies the name of the WVAR objects.
#' @param cex_labels A \code{double} that specifies the magnification of the labels (x and y).
#' @author Stephane Guerrier and Justin Lee
#' @export
#' @examples
#' set.seed(999)
#' n = 10^4
#' Xt = arima.sim(n = n, list(ar = 0.10))
#' Yt = arima.sim(n = n, list(ar = 0.35))
#' Zt = arima.sim(n = n, list(ar = 0.70))
#' Wt = arima.sim(n = n, list(ar = 0.95))
#'
#' wv_Xt = wvar(Xt)
#' wv_Yt = wvar(Yt)
#' wv_Zt = wvar(Zt)
#' wv_Wt = wvar(Wt)
#'
#' compare_wvar(wv_Xt, wv_Yt, wv_Zt, wv_Wt)
compare_wvar = function(... , split = FALSE, add_legend = TRUE, units = NULL, xlab = NULL,
ylab = NULL, main = NULL, col_wv = NULL, col_ci = NULL, nb_ticks_x = NULL,
nb_ticks_y = NULL, legend_position = NULL, ci_wv = NULL, point_cex = NULL,
point_pch = NULL, names = NULL, cex_labels = 0.8){
obj_list = list(...)
obj_name = as.character(substitute(...()))
obj_len = length(obj_list)
# Check if passed objects are of class wvar
is_wvar = sapply(obj_list, FUN = is, class2 = 'wvar')
if(!all(is_wvar == T)){
stop("Supplied objects must be 'wvar' objects.")
}
# Check length of time series argument
if (obj_len == 0){
stop('No object given!')
}else if (obj_len == 1){
# -> plot.wvar
plot.wvar(..., nb_ticks_x = nb_ticks_x, nb_ticks_y = nb_ticks_y)
}else{
if (is.null(xlab)){
if (is.null(units)){
xlab = expression(paste("Scale ", tau, sep =""))
}else{
xlab = bquote(paste("Scale ", "(", .(units), ")", sep = " "))
}
}else{
xlab = xlab
}
if (is.null(ylab)){
ylab = bquote(paste("Wavelet Variance ", nu^2, sep = " "))
}else{
ylab = ylab
}
if (is.null(ci_wv)){
ci_wv = rep(TRUE, obj_len)
}else{
ci_wv = rep(ci_wv, obj_len)
}
# Main Title
if (split == FALSE){
if (is.null(main)){
main = "Haar Wavelet Variance Representation"
}
}else{
if (!is.null(main) && (dim(main)[1] != obj_len || dim(main)[2] != obj_len)){
main = NULL
}
}
hues = seq(15, 375, length = obj_len + 1)
# Line Colors
if (is.null(col_wv)){
col_wv = hcl(h = hues, l = 65, c = 200, alpha = 1)[seq_len(obj_len)]
}else{
if (length(col_wv) != obj_len){
col_wv = hcl(h = hues, l = 65, c = 200, alpha = 1)[seq_len(obj_len)]
}
}
# CI Colors
if (is.null(col_ci)){
col_ci = hcl(h = hues, l = 80, c = 100, alpha = 0.2)[seq_len(obj_len)]
}else{
if (length(col_ci) != obj_len){
col_ci = hcl(h = hues, l = 80, c = 100, alpha = 0.2)[seq_len(obj_len)]
}
}
# X and Y Limits
x_range = y_range = rep(NULL, 2)
for (i in 1:obj_len){
x_range = range(c(x_range, obj_list[[i]]$scales))
y_range = range(c(y_range, obj_list[[i]]$ci_low, obj_list[[i]]$ci_high))
}
x_low = floor(log10(x_range[1]))
x_high = ceiling(log10(x_range[2]))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
# Axes Labels and Ticks
if (is.null(nb_ticks_x)){
nb_ticks_x = 6
}
if (is.null(nb_ticks_y)){
nb_ticks_y = 5
}
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > nb_ticks_x){
x_ticks = x_low + ceiling((x_high - x_low)/(nb_ticks_x + 1))*(0:nb_ticks_x)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(10^ .(i))))
x_at = 10^x_ticks
x_actual_length = sum((x_at < x_range[2])*(x_at > x_range[1]))
if (x_actual_length < (3 + as.numeric(split == FALSE))){
x_low = floor(log2(x_range[1]))
x_high = ceiling(log2(x_range[2]))
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > 8){
x_ticks = seq(x_low, x_high, by = 2)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(2^ .(i))))
x_at = 2^x_ticks
}
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels = sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
y_at = 10^y_ticks
# Legend position
if (is.null(legend_position)){
inter = rep(NA, obj_len)
for (i in 1:obj_len){
inter[i] = obj_list[[i]]$variance[1]
}
mean_wv_1 = mean(inter)
if (which.min(abs(c(y_low, y_high) - log2(mean_wv_1))) == 1){
legend_position = "topleft"
}else{
legend_position = "bottomleft"
}
}
# Type of Points
if (is.null(point_pch)){
inter = rep(15:18, obj_len)
point_pch = inter[1:obj_len]
}else{
if (length(point_pch) != obj_len){
inter = rep(15:18, obj_len)
point_pch = inter[1:obj_len]
}
}
#Size of Points
if (is.null(point_cex)){
inter = rep(c(1.25,1.25,1.25,1.25), obj_len)
point_cex = inter[1:obj_len]
}else{
if (length(point_pch) != obj_len){
inter = rep(c(1.25,1.25,1.25,1.25), obj_len)
point_cex = inter[1:obj_len]
}
}
# Names of WVAR Objects
if (is.null(names)){
names = obj_name
}else{
if (length(names) != obj_len){
names = obj_name
}
}
# Arguments passed into compare_wvar_split or compare_wvar_no_split
graph_details = list(obj_list = obj_list, obj_len = obj_len, names = names, xlab = xlab,
ylab = ylab, col_wv = col_wv, add_legend = add_legend,
col_ci = col_ci, main = main, legend_position = legend_position,
ci_wv = ci_wv, point_cex = point_cex, point_pch = point_pch,
x_range = x_range, y_range = y_range, x_ticks = x_ticks,
x_labels = x_labels, y_labels = y_labels, x_at = x_at, y_at = y_at,
y_ticks = y_ticks, nb_ticks_x = nb_ticks_x, nb_ticks_y = nb_ticks_y,
cex_labels = cex_labels)
if (split == FALSE){
# -> compare_wvar_no_split
compare_wvar_no_split(graph_details)
}else{
# -> compare_wvar_split
compare_wvar_split(graph_details)
}
}
}
Try the wv package in your browser
Any scripts or data that you put into this service are public.
wv documentation built on Jan. 17, 2020, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions unitConversion wvar wvar.lts wvar.gts wvar.ts wvar.imu wvar.default create_wvar print.wvar summary.wvar plot.wvar
Documented in create_wvar plot.wvar print.wvar summary.wvar unitConversion wvar wvar.default wvar.gts wvar.imu wvar.lts wvar.ts
#' @title Convert Unit of Time Series Data
#' @description Manipulate the units of time to different ones
#' @keywords internal
#' @param x A \code{vector} containing the values on x-axis.
#' @param from.unit A \code{string} indicating the unit which the data is converted from.
#' @param to.unit A \code{string} indicating the unit which the data is converted to.
#' @details
#' The supported units are "ns"(nanosecond), "ms"(millisecond), "sec", "min", "hour", "day", "month", and "year".
#' Make sure \code{from.unit} and \code{to.unit} are not \code{NULL} before it is passed to this function.
#' @return A \code{list} with the following structure:
#' \itemize{
#' \item "x": Data
#' \item "converted": A \code{boolean} indicating whether conversion is made
#' }
#' @export
#' @examples
#' x = seq(60, 3600, 60)
#' unitConversion(x, 'sec', 'min')
#' y = 1:10
#' unitConversion(y, 'hour', 'sec')
unitConversion = function(x, from.unit, to.unit){
#ns, ms, second, min, hour, day, month, year
unit = c(ns = 1, ms = 2,se = 3, mi = 4, ho = 5, da = 6, mo = 7, ye = 8)
#assume 1 month = 30 days
ratio = c(1E6, 1E3, 60, 60, 24, 30, 12)
from.unit.1 = substr(from.unit, 1, 2)
to.unit.1 = substr(to.unit, 1, 2)
#check unit:
no.convert = F
if(from.unit.1 == to.unit.1){no.convert = T}
if(is.na(unit[from.unit.1]) ) {
message = paste('No such unit: ', from.unit, '. Supported units are "ns"(nanosecond), "ms"(millisecond), "sec", "min", "hour", "day", "month", and "year". Conversion is terminated.', sep = '')
warning(message); no.convert = T}
if(is.na(unit[to.unit.1]) ) {
message = paste('No such unit: ', to.unit, '. Supported units are "ns"(nanosecond), "ms"(millisecond), "sec", "min", "hour", "day", "month", and "year". Conversion is terminated.', sep = '')
warning(message); no.convert = T}
if(!no.convert){
#print out warning when day is convert to month, or month is converted to day.
conversionRange = unit[from.unit.1] : unit[to.unit.1]
if(6 %in% conversionRange && 7 %in% conversionRange){
warning('Unit conversion might be wrong because this function simply assumes 1 month = 30 days.')
}
}
if(!no.convert){
if(unit[from.unit.1] > unit[to.unit.1]){
temp = ratio[unit[to.unit.1]: (unit[from.unit.1]-1)]
multiplier = prod(temp)
x = x*multiplier
}else{
temp = ratio[unit[from.unit.1]: (unit[to.unit.1]-1) ]
multiplier = prod(temp)
x = x/multiplier
}
}
obj = list(x = x, converted = !no.convert)
return(obj)
}
#' @title Wavelet Variance
#' @description Calculates the (MO)DWT wavelet variance
#' @param x A \code{vector} with dimensions N x 1.
#' @param decomp A \code{string} that indicates whether to use a "dwt" or "modwt" decomposition.
#' @param filter A \code{string} that specifies which wavelet filter to use.
#' @param nlevels An \code{integer} that indicates the level of decomposition. It must be less than or equal to floor(log2(length(x))).
#' @param alpha A \code{double} that specifies the significance level which in turn specifies the \eqn{1-\alpha} confidence level.
#' @param robust A \code{boolean} that triggers the use of the robust estimate.
#' @param eff A \code{double} that indicates the efficiency as it relates to an MLE.
#' @param freq A \code{numeric} that provides the rate of samples.
#' @param from.unit A \code{string} indicating the unit from which the data is converted.
#' @param to.unit A \code{string} indicating the unit to which the data is converted.
#' @param ... Further arguments passed to or from other methods.
#' @return A \code{list} with the structure:
#' \itemize{
#' \item "variance": Wavelet Variance
#' \item "ci_low": Lower CI
#' \item "ci_high": Upper CI
#' \item "robust": Robust active
#' \item "eff": Efficiency level for Robust calculation
#' \item "alpha": p value used for CI
#' \item "unit": String representation of the unit
#' }
#' @details
#' The default value of \code{nlevels} will be set to \eqn{\left\lfloor {{{\log }_2}\left( {length\left( x \right)} \right)} \right\rfloor}{floor(log2(length(x)))}, unless otherwise specified.
#' @author James Balamuta, Justin Lee and Stephane Guerrier
#' @rdname wvar
#' @examples
#' set.seed(999)
#' x = rnorm(100)
#'
#' # Default
#' wvar(x)
#'
#' # Robust
#' wvar(x, robust = TRUE, eff=0.3)
#'
#' # Classical
#' wvar(x, robust = FALSE, eff=0.3)
#'
#' # 90% Confidence Interval
#' wvar(x, alpha = 0.10)
#' @export
wvar = function(x, ...) {
if (sum(class(x) %in% "gts") == 1){
x = as.numeric(x)
}
UseMethod("wvar")
}
#' @rdname wvar
#' @export
wvar.lts = function(x, decomp = "modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, to.unit = NULL, ...){
warning('`lts` object is detected. This function can only operate on the combined process.')
freq = attr(x, 'freq')
unit = attr(x, 'unit')
x = x[,ncol(x)]
wvar.default(x, decomp, filter, nlevels, alpha, robust, eff, freq = freq, from.unit = unit, to.unit = to.unit)
}
#' @rdname wvar
#' @export
wvar.gts = function(x, decomp="modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, to.unit = NULL, ...){
freq = attr(x, 'freq')
unit = attr(x, 'unit')
x = x[,1]
wvar.default(x, decomp, filter, nlevels, alpha, robust, eff, freq = freq, from.unit = unit, to.unit = to.unit)
}
#' @rdname wvar
#' @export
wvar.ts = function(x, decomp="modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, to.unit = NULL, ...){
freq = attr(x, 'tsp')[3]
unit = NULL
wvar.default(x, decomp, filter, nlevels, alpha, robust, eff, freq = freq, from.unit = unit, to.unit = to.unit)
}
#' @rdname wvar
#' @export
wvar.imu = function(x, decomp="modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, to.unit = NULL, ...){
# Retrive sensor name
if (!is.null(attr(x, "stype"))){
sensor_name = attr(x, "stype")
}else{
warning("Unknown sensor name. IMU object is missing some information.")
sensor_name = NULL
}
# Retrive freq
if (!is.null(attr(x, "freq"))){
freq = attr(x, "freq")
}else{
warning("Unknown frequency. IMU object is missing some information. Freq is set to 1 by default.")
freq = 1
}
# Retrive sample size
if (!is.null(attr(x, "dim"))){
n = attr(x, "dim")[1]
}else{
warning("Unknown sample size. IMU object is missing some information.")
n = NULL
}
# Retrive col names
if (!is.null(attr(x, "dimnames")[[2]])){
col_names = attr(x, "dimnames")[[2]]
}else{
stop("Unknown colunms names. IMU object is missing some information.")
col_names = NULL
}
# Retrive sensor
if (!is.null(attr(x, "sensor"))){
sensor = attr(x, "sensor")
}else{
warning("Unknown sensor. IMU object is missing some information.")
sensor = NULL
}
# Retrive axis
if (!is.null(attr(x, "axis"))){
ax = attr(x, "axis")
}else{
warning("Unknown axes. IMU object is missing some information.")
ax = NULL
}
# Compute wvar
m = length(col_names)
wvariance = list()
for (i in 1:m){
wvariance[[i]] = wvar.default(x[,i], decomp, filter, nlevels, alpha, robust, eff, freq = freq, to.unit = to.unit)
}
names(wvariance) = col_names
out = list(sensor = sensor_name, freq = freq, n = n, type = sensor, axis = ax, wvar = wvariance)
class(out) = "imu_wvar"
invisible(out)
}
#' @rdname wvar
#' @export
#' @importFrom methods is
wvar.default = function(x, decomp = "modwt", filter = "haar", nlevels = NULL, alpha = 0.05, robust = FALSE, eff = 0.6, freq = 1, from.unit = NULL, to.unit = NULL, ...){
if(is.null(x)){
stop("`x` must contain a value")
}else if((is.data.frame(x) || is.matrix(x))){
if(ncol(x) > 1) stop("There must be only one column of data supplied.")
}
if(decomp == "modwt" && is.null(nlevels)){
nlevels = floor(log2(length(x))-1)
}else if(decomp == "dwt" && is.null(nlevels)){
nlevels = floor(log2(length(x)))
}
# Check Freq
if(!is(freq,"numeric") || length(freq) != 1){ stop("'freq' must be one numeric number.") }
if(freq <= 0) { stop("'freq' must be larger than 0.") }
# Check Unit
all.units = c('ns', 'ms', 'sec', 'second', 'min', 'minute', 'hour', 'day', 'mon', 'month', 'year')
if( (!is.null(from.unit) && !from.unit %in% all.units) || (!is.null(to.unit) && !to.unit %in% all.units) ){
stop('The supported units are "ns", "ms", "sec", "min", "hour", "day", "month", "year". ')
}
if(robust) {
if(eff > 0.99) {
stop("The efficiency specified is too close to the classical case. Use `robust = FALSE`")
}
}
obj = modwt_wvar_cpp(signal=x, nlevels=nlevels, robust=robust, eff=eff, alpha=alpha,
ci_type="eta3", strWavelet=filter, decomp = decomp)
# nlevels may be changed during modwt
nlevels = nrow(obj)
scales = scales_cpp(nlevels)/freq
# NO Unit Conversion
if( is.null(from.unit) && is.null(to.unit)==F ){
warning("'from.unit' is NULL. Unit conversion was not done.")
}
# Unit Conversion
if (!is.null(from.unit)){
if (!is.null(to.unit)){
convert.obj = unitConversion(scales, from.unit = from.unit, to.unit = to.unit)
if (convert.obj$converted) {
# YES Unit Conversion
scales = convert.obj$x
message(paste0('Unit of object is converted from ', from.unit, ' to ', to.unit), appendLF = T)
}
}
}
if(!is.null(from.unit) && !is.null(to.unit)){
unit = to.unit
}else{
unit = from.unit}
create_wvar(obj, decomp, filter, robust, eff, alpha, scales, unit)
}
#' @title Create a \code{wvar} object
#' @description Structures elements into a \code{wvar} object
#' @param obj A \code{matrix} with dimensions N x 3 that contains Wavelet Variance, Lower CI, and Upper CI.
#' @param decomp A \code{string} that indicates whether to use a "dwt" or "modwt" decomposition.
#' @param filter A \code{string} that specifies the type of wavelet filter used in the decomposition.
#' @param robust A \code{boolean} that triggers the use of the robust estimate.
#' @param eff A \code{double} that indicates the efficiency as it relates to an MLE.
#' @param alpha A \code{double} that specifies the significance level which in turn specifies the \eqn{1-\alpha} confidence level.
#' @param scales A \code{vec} that contains the amount of decomposition performed at each level.
#' @param unit A \code{string} that indicates the unit expression of the frequency.
#' @return A \code{list} with the structure:
#' \itemize{
#' \item "variance": Wavelet variance
#' \item "ci_low": Lower CI
#' \item "ci_high": Upper CI
#' \item "robust": Robust active
#' \item "eff": Efficiency level for robust calculation
#' \item "alpha": p value used for CI
#' \item "unit": String representation of the unit
#' }
#' @keywords internal
create_wvar = function(obj, decomp, filter, robust, eff, alpha, scales, unit){
structure(list(variance = obj[,1],
ci_low = obj[,2],
ci_high = obj[,3],
robust = robust,
eff = eff,
alpha = alpha,
scales = scales,
decomp = decomp,
unit = unit,
filter = filter), class = "wvar")
}
#' @title Print Wavelet Variances
#' @description Displays the summary table of wavelet variance.
#' @author James Balamuta
#' @method print wvar
#' @export
#' @keywords internal
#' @param x A \code{wvar} object.
#' @param ... Further arguments passed to or from other methods.
#' @return Summary table
#' @examples
#' set.seed(999)
#' x = rnorm(100)
#' out = wvar(x)
#' print( out )
print.wvar = function(x, ...){
mat = matrix(unlist(x[1:3]),ncol=3,byrow=F)
colnames(mat) = c("Variance", "Low CI", "High CI")
rownames(mat) = x$scales
print(mat)
}
#' @title Summary of Wavelet Variances
#' @description Displays the summary table of wavelet variance accounting for CI values and supplied efficiency.
#' @method summary wvar
#' @export
#' @keywords internal
#' @param object A \code{wvar} object.
#' @param ... Additional arguments affecting the summary produced.
#' @return Summary table and other properties of the object.
#' @author James Balamuta
#' @examples
#' set.seed(999)
#' x = rnorm(100)
#' ret = wvar(x)
#' summary(ret)
summary.wvar = function(object, ...){
name = if(object$robust){
"robust"
}else{
"classical"
}
cat("Results of the wavelet variance calculation using the ",name, " method.\n",sep="")
if(object$robust){
cat("Robust was created using efficiency=",object$eff,"\n",sep="")
}
cat("The confidence interval was generated using (1-",object$alpha,")*100 \n",sep="")
print(object)
}
#' @title Plot Wavelet Variance
#' @description Displays a plot of wavelet variance accounting for CI values and supplied efficiency.
#' @method plot wvar
#' @keywords internal
#' @param x A \code{wvar} object.
#' @param units A \code{string} that specifies the units of time plotted on the x axis.
#' @param xlab A \code{string} that gives a title for the x axis.
#' @param ylab A \code{string} that gives a title for the y axis.
#' @param main A \code{string} that gives an overall title for the plot.
#' @param col_wv A \code{string} that specifies the color of the wavelet variance line.
#' @param col_ci A \code{string} that specifies the color of the confidence interval polygon.
#' @param nb_ticks_x An \code{integer} that specifies the maximum number of ticks for the x-axis.
#' @param nb_ticks_y An \code{integer} that specifies the maximum number of ticks for the y-axis.
#' @param legend_position A \code{string} that specifies the position of the legend (use \code{legend_position = NA} to remove legend).
#' @param ci_wv A \code{boolean} that determines whether a confidence interval polygon will be drawn.
#' @param point_cex A \code{double} that specifies the size of each symbol to be plotted.
#' @param point_pch A \code{double} that specifies the symbol type to be plotted.
#' @param ... Additional arguments affecting the plot.
#' @return Plot of wavelet variance and confidence interval for each scale.
#' @author Stephane Guerrier, Nathanael Claussen, and Justin Lee
#' @export
#' @examples
#' set.seed(999)
#' n = 10^4
#' Xt = rnorm(n)
#' wv = wvar(Xt)
#' plot(wv)
#' plot(wv, main = "Simulated white noise", xlab = "Scales")
#' plot(wv, units = "sec", legend_position = "topright")
#' plot(wv, col_wv = "darkred", col_ci = "pink")
plot.wvar = function(x, units = NULL, xlab = NULL, ylab = NULL, main = NULL,
col_wv = NULL, col_ci = NULL, nb_ticks_x = NULL, nb_ticks_y = NULL,
legend_position = NULL, ci_wv = NULL, point_cex = NULL,
point_pch = NULL, ...){
# Labels
if (is.null(xlab)){
if (is.null(units)){
xlab = expression(paste("Scale ", tau, sep =""))
}else{
xlab = bquote(paste("Scale ", tau, " [", .(units), "]", sep = " "))
}
}
if (is.null(ylab)){
ylab = expression(paste("Wavelet Variance ", nu^2, sep = ""))
}else{
ylab = ylab
}
# Main Title
if (is.null(main)){
main = "Haar Wavelet Variance Representation"
}
# Line and CI colors
if (is.null(col_wv)){
col_wv = "darkblue"
}
if (is.null(col_ci)){
col_ci = hcl(h = 210, l = 65, c = 100, alpha = 0.2)
}
# Range
x_range = range(x$scales)
x_low = floor(log2(x_range[1]))
x_high = ceiling(log2(x_range[2]))
y_range = range(c(x$ci_low, x$ci_high))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
# Axes
if (is.null(nb_ticks_x)){
nb_ticks_x = 6
}
if (is.null(nb_ticks_y)){
nb_ticks_y = 5
}
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > nb_ticks_x){
x_ticks = x_low + ceiling((x_high - x_low)/(nb_ticks_x + 1))*(0:nb_ticks_x)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(2^ .(i))))
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels <- sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
# Legend Position
if (is.null(legend_position)){
if (which.min(abs(c(y_low, y_high) - log2(x$variance[1]))) == 1){
legend_position = "topleft"
}else{
legend_position = "bottomleft"
}
}
# Main Plot
plot(NA, xlim = x_range, ylim = y_range, xlab = xlab, ylab = ylab,
log = "xy", xaxt = 'n', yaxt = 'n', bty = "n", ann = FALSE)
win_dim = par("usr")
par(new = TRUE)
plot(NA, xlim = x_range, ylim = 10^c(win_dim[3], win_dim[4] + 0.09*(win_dim[4] - win_dim[3])),
xlab = xlab, ylab = ylab, log = "xy", xaxt = 'n', yaxt = 'n', bty = "n")
win_dim = par("usr")
# Add Grid
abline(v = 2^x_ticks, lty = 1, col = "grey95")
abline(h = 10^y_ticks, lty = 1, col = "grey95")
# Add Title
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
polygon(x_vec, y_vec, col = "grey95", border = NA)
text(x = 10^mean(c(win_dim[1], win_dim[2])), y = 10^(win_dim[4] - 0.09/2*(win_dim[4] - win_dim[3])), main)
# Add Axes and Box
lines(x_vec[1:2], rep(10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])),2), col = 1)
#y_ticks = y_ticks[(2^y_ticks) < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))]
y_labels = y_labels[1:length(y_ticks)]
box()
axis(1, at = 2^x_ticks, labels = x_labels, padj = 0.3)
axis(2, at = 10^y_ticks, labels = y_labels, padj = -0.2)
# CI for WV
if (ci_wv == TRUE || is.null(ci_wv)){
polygon(c(x$scales, rev(x$scales)), c(x$ci_low, rev(x$ci_high)),
border = NA, col = col_ci)
}
# Add legend
CI_conf = 1 - x$alpha
if (x$robust == TRUE){
wv_title_part1 = "Empirical Robust WV "
}else{
wv_title_part1 = "Empirical WV "
}
if (!is.na(legend_position)){
if (legend_position == "topleft"){
legend_position = 10^c(1.1*win_dim[1], 0.98*(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])))
legend(x = legend_position[1], y = legend_position[2],
legend = c(as.expression(bquote(paste(.(wv_title_part1), hat(nu)^2))),
as.expression(bquote(paste("CI(",hat(nu)^2,", ",.(CI_conf),")")))),
pch = c(16, 15), lty = c(1, NA), col = c(col_wv, col_ci), cex = 1, pt.cex = c(1.25, 3), bty = "n")
}else{
if (legend_position == "topright"){
legend_position = 10^c(0.7*win_dim[2], 0.98*(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])))
legend(x = legend_position[1], y = legend_position[2],
legend = c(as.expression(bquote(paste(.(wv_title_part1), hat(nu)^2))),
as.expression(bquote(paste("CI(",hat(nu)^2,", ",.(CI_conf),")")))),
pch = c(16, 15), lty = c(1, NA), col = c(col_wv, col_ci), cex = 1, pt.cex = c(1.25, 3), bty = "n")
}else{
legend(legend_position,
legend = c(as.expression(bquote(paste(.(wv_title_part1), hat(nu)^2))),
as.expression(bquote(paste("CI(",hat(nu)^2,", ",.(CI_conf),")")))),
pch = c(16, 15), lty = c(1, NA), col = c(col_wv, col_ci), cex = 1, pt.cex = c(1.25, 3), bty = "n")
}
}
}
# Add WV
lines(x$scales, x$variance, type = "l", col = col_wv, pch = 16)
if (is.null(point_pch)){
point_pch = 16
}
if (is.null(point_cex)){
point_cex = 1.25
}
lines(x$scales, x$variance, type = "p", col = col_wv, pch = point_pch, cex = point_cex)
}
#' @title Plot Wavelet Variance based on IMU Data
#' @description Displays a plot of wavelet variance accounting for CI values and supplied efficiency.
#' @method plot imu_wvar
#' @keywords internal
#' @param x A \code{wvar} object.
#' @param xlab A \code{string} that gives a title for the x axis.
#' @param ylab A \code{string} that gives a title for the y axis.
#' @param main A \code{string} that gives an overall title for the plot.
#' @param col_wv A \code{string} that specifies the color of the wavelet variance line.
#' @param col_ci A \code{string} that specifies the color of the confidence interval polygon.
#' @param nb_ticks_x An \code{integer} that specifies the maximum number of ticks for the x-axis.
#' @param nb_ticks_y An \code{integer} that specifies the maximum number of ticks for the y-axis.
#' @param ci_wv A \code{boolean} that determines whether a confidence interval polygon will be drawn.
#' @param point_cex A \code{double} that specifies the size of each symbol to be plotted.
#' @param point_pch A \code{double} that specifies the symbol type to be plotted.
#' @param ... Additional arguments affecting the plot.
#' @return Plot of wavelet variance and confidence interval for each scale.
#' @author Stephane Guerrier and Yuming Zhang
#' @export
#' @examples
#' data("kvh1750_wv")
#' plot(kvh1750_wv)
plot.imu_wvar = function(x, xlab = NULL, ylab = NULL, main = NULL,
col_wv = NULL, col_ci = NULL, nb_ticks_x = NULL, nb_ticks_y = NULL,
ci_wv = NULL, point_cex = NULL, point_pch = NULL, ...){
type = unique(x$type)
if ("Gyroscope" %in% type){
gyro_index = which(x$type == "Gyroscope")
}else{
gyro_index = NULL
}
if ("Accelerometer" %in% type){
accel_index = which(x$type == "Accelerometer")
}else{
accel_index = NULL
}
ncol = length(unique(x$axis))
nrow = length(type)
m = length(x$wvar)
J = length(x$wvar[[1]]$variance)
# remove negative CI values
index_to_remove = c()
for (i in 1:m) {
if(length(which(x$wvar[[i]]$lci<0)) > 0){
index_to_remove = c(index_to_remove, which(x$wvar[[i]]$lci<0))
}
}
if (!is.null(index_to_remove)){
index_to_remove = unique(index_to_remove)
index_to_keep = which(seq(1:J) != index_to_remove)
}else{
index_to_keep = 1:J
}
J = length(index_to_keep)
scales = x$wvar[[1]]$scales[index_to_keep]
ci_up = ci_lw = av = matrix(NA, J, m)
for (i in 1:m){
ci_up[,i] = x$wvar[[i]]$ci_high[index_to_keep]
ci_lw[,i] = x$wvar[[i]]$ci_low[index_to_keep]
av[,i] = x$wvar[[i]]$variance[index_to_keep]
}
# Axes
if (is.null(nb_ticks_x)){
nb_ticks_x = 6
}
if (is.null(nb_ticks_y)){
nb_ticks_y = 5
}
# Range
x_range = range(scales)
x_low = floor(log10(x_range[1]))
x_high = ceiling(log10(x_range[2]))
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > nb_ticks_x){
x_ticks = x_low + ceiling((x_high - x_low)/(nb_ticks_x + 1))*(0:nb_ticks_x)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(10^ .(i))))
# Line and CI colors
if (is.null(col_wv)){
col_wv = "darkblue"
}
if (is.null(col_ci)){
col_ci = hcl(h = 210, l = 65, c = 100, alpha = 0.2)
}
if (is.null(point_pch)){
point_pch = 16
}
if (is.null(point_cex)){
point_cex = 1.25
}
# Main Title
if (is.null(main)){
main = paste("Wavelet Variance Representation - ", x$sensor, " @ ", x$freq, " Hz", sep="")
}
# Labels
if (is.null(xlab)){
xlab = bquote(paste("Averaging time ", tau, " [sec]", sep = " "))
}
if (is.null(ylab)){
ylab = expression(paste("Wavelet Variance ", nu, sep = ""))
}
# Main plot
par(omi=rep(1, 4), mar=c(0,0,0,0), mfrow=c(nrow,ncol))
# Gyro
if (!is.null(gyro_index)){
y_range = c(min(ci_lw[,gyro_index]), max(ci_up[,gyro_index]))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels <- sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
for (i in seq_along(gyro_index)){
plot(NA, xlim = range(scales), ylim = y_range, xaxt="n", yaxt="n", log = "xy", bty = "n")
box(col = "grey")
mtext(paste("Axis - ", x$axis[gyro_index][i], sep = ""), 3, line = 0.5)
if (i == 1){
axis(2, at = 10^y_ticks, labels = y_labels, padj = -0.2, cex = 1.25)
}
if (i == 1){
mtext("Gyroscope", 2, line = 4.5)
mtext(ylab, 2, line = 2.5)
}
abline(h = 10^y_ticks, col = "grey85")
abline(v = 10^x_ticks, col = "grey85")
# CI for AD
if(ci_wv == TRUE || is.null(ci_wv)){
polygon(c(scales, rev(scales)), c(ci_lw[,gyro_index[i]], rev(ci_up[,gyro_index[i]])),
border = NA, col = col_ci)
}
# Add AD
lines(scales, (av[,gyro_index[i]]), type = "l", col = col_wv, pch = 16)
lines(scales, (av[,gyro_index[i]]), type = "p", col = col_wv, pch = point_pch, cex = point_cex)
if (is.null(accel_index)){
axis(1, at = 10^x_ticks, labels = x_labels, padj = -0.2, cex = 1.25)
}
}
}
# Accel
if (!is.null(accel_index)){
y_range = c(min(ci_lw[,accel_index]), max(ci_up[,accel_index]))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels <- sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
for (i in seq_along(accel_index)){
plot(NA, xlim = range(scales), ylim = y_range, xaxt="n", yaxt="n", log = "xy", bty = "n")
box(col = "grey")
if (i == 1){
axis(2, at = 10^y_ticks, labels = y_labels, padj = -0.2, cex = 1.25)
}
if (i == 1){
mtext("Accelerometer", 2, line = 4.5)
mtext(ylab, 2, line = 2.5)
}
if (length(accel_index) == 3 && i == 2){
mtext(xlab, 1, line = 3.5)
}
if (is.null(gyro_index)){
mtext(paste("Axis - ", x$axis[gyro_index][i], sep = ""), 3, line = 0.5)
}
abline(h = 10^y_ticks, col = "grey85")
abline(v = 10^x_ticks, col = "grey85")
# CI for AD
if(ci_wv == TRUE || is.null(ci_wv)){
polygon(c(scales, rev(scales)), c(ci_lw[,accel_index[i]], rev(ci_up[,accel_index[i]])),
border = NA, col = col_ci)
}
# Add AD
lines(scales, (av[,accel_index[i]]), type = "l", col = col_wv, pch = 16)
lines(scales, (av[,accel_index[i]]), type = "p", col = col_wv, pch = point_pch, cex = point_cex)
axis(1, at = 10^x_ticks, labels = x_labels, padj = -0.2, cex = 1.25)
}
}
# Add main title
mtext(main, side = 3, line = 3, outer = TRUE)
par(mfrow = c(1,1))
}
#' @title Comparison between classical and robust Wavelet Variances
#' @description Displays a plot of the wavelet variances (classical and robust) for a given time series accounting for CI values.
#' @param x A time series objects.
#' @param eff An \code{integer} that specifies the efficiency of the robust estimator.
#' @param units A \code{string} that specifies the units of time plotted on the x axis.
#' @param xlab A \code{string} that gives a title for the x axis.
#' @param ylab A \code{string} that gives a title for the y axis.
#' @param main A \code{string} that gives an overall title for the plot.
#' @param col_wv A \code{string} that specifies the color of the wavelet variance line.
#' @param col_ci A \code{string} that specifies the color of the confidence interval shade.
#' @param nb_ticks_x An \code{integer} that specifies the maximum number of ticks for the x-axis.
#' @param nb_ticks_y An \code{integer} that specifies the maximum number of ticks for the y-axis.
#' @param legend_position A \code{string} that specifies the position of the legend (use \code{legend_position = NA} to remove legend).
#' @param ... Additional arguments affecting the plot.
#' @return Plot of wavelet variance and confidence interval for each scale.
#' @author Stephane Guerrier, Nathanael Claussen, and Justin Lee
#' @examples
#' set.seed(999)
#' n = 10^4
#' Xt = rnorm(n)
#' wv = wvar(Xt)
#'
#' plot(wv)
#' plot(wv, main = "Simulated white noise", xlab = "Scales")
#' plot(wv, units = "sec", legend_position = "topright")
#' plot(wv, col_wv = "darkred", col_ci = "pink")
#' @export
robust_eda = function(x, eff = 0.6, units = NULL, xlab = NULL, ylab = NULL, main = NULL,
col_wv = NULL, col_ci = NULL, nb_ticks_x = NULL, nb_ticks_y = NULL,
legend_position = NULL, ...){
wv_cl = wvar(x)
wv_rob = wvar(x, robust = TRUE, eff = eff)
# Labels
if (is.null(xlab)){
if (is.null(units)){
xlab = expression(paste("Scale ", tau, sep =""))
}else{
xlab = bquote(paste("Scale ", "", tau , " [", .(units), "]", sep = ""))
}
}
if (is.null(ylab)){
ylab = expression(paste("Wavelet Variance ", nu^2, sep = ""))
}else{
ylab = ylab
}
# Main Title
if (is.null(main)){
main = "Classical vs Robust WV"
}
# Line and CI colors
if (is.null(col_wv)){
col_wv = c("darkblue", "darkorange2")
}
if (is.null(col_ci)){
col_ci = c(hcl(h = 210, l = 65, c = 100, alpha = 0.2), hcl(h = 60, l = 65, c = 100, alpha = 0.2))
}
# Range
x_range = range(wv_cl$scales)
x_low = floor(log2(x_range[1]))
x_high = ceiling(log2(x_range[2]))
y_range = range(c(wv_cl$ci_low, wv_cl$ci_high, wv_rob$ci_low, wv_rob$ci_high))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
# Axes
if (is.null(nb_ticks_x)){
nb_ticks_x = 6
}
if (is.null(nb_ticks_y)){
nb_ticks_y = 5
}
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > nb_ticks_x){
x_ticks = x_low + ceiling((x_high - x_low)/(nb_ticks_x + 1))*(0:nb_ticks_x)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(2^ .(i))))
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels <- sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
# Legend position
if (is.null(legend_position)){
if (which.min(abs(c(y_low, y_high) - log2(wv_rob$variance[1]))) == 1){
legend_position = "topleft"
}else{
legend_position = "bottomleft"
}
}
# Main plot
plot(NA, xlim = x_range, ylim = y_range, xlab = xlab, ylab = ylab,
log = "xy", xaxt = 'n', yaxt = 'n', bty = "n", ann = FALSE)
win_dim = par("usr")
par(new = TRUE)
plot(NA, xlim = x_range, ylim = 10^c(win_dim[3], win_dim[4] + 0.09*(win_dim[4] - win_dim[3])),
xlab = xlab, ylab = ylab, log = "xy", xaxt = 'n', yaxt = 'n', bty = "n")
win_dim = par("usr")
# Add grid
abline(v = 2^x_ticks, lty = 1, col = "grey95")
abline(h = 10^y_ticks, lty = 1, col = "grey95")
# Add title
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
polygon(x_vec, y_vec, col = "grey95", border = NA)
text(x = 10^mean(c(win_dim[1], win_dim[2])), y = 10^(win_dim[4] - 0.09/2*(win_dim[4] - win_dim[3])), main)
# Add Axes and Box
lines(x_vec[1:2], rep(10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])),2), col = 1)
y_ticks = y_ticks[(2^y_ticks) < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))]
y_labels = y_labels[1:length(y_ticks)]
box()
axis(1, at = 2^x_ticks, labels = x_labels, padj = 0.3)
axis(2, at = 10^y_ticks, labels = y_labels, padj = -0.2)
# CI for WV
polygon(c(wv_cl$scales, rev(wv_cl$scales)), c(wv_cl$ci_low, rev(wv_cl$ci_high)),
border = NA, col = col_ci[1])
polygon(c(wv_rob$scales, rev(wv_rob$scales)), c(wv_rob$ci_low, rev(wv_rob$ci_high)),
border = NA, col = col_ci[2])
# Legend Position
if (!is.na(legend_position)){
if (legend_position == "topleft"){
legend_position = 10^c(1.1*win_dim[1], 0.98*(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])))
legend(x = legend_position[1], y = legend_position[2],
legend = c("Classical WV", "Classical CI", "Robust WV", "Robust CI"),
pch = c(16, 15, 16, 15), lty = c(1, NA, 1, NA), col = c(col_wv[1], col_ci[1], col_wv[2], col_ci[2]),
cex = 1, pt.cex = c(1.25, 3, 1.25, 3), bty = "n")
}else{
if (legend_position == "topright"){
legend_position = 10^c(0.7*win_dim[2], 0.98*(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])))
legend(x = legend_position[1], y = legend_position[2],
legend = c("Classical WV", "Classical CI", "Robust WV", "Robust CI"),
pch = c(16, 15, 16, 15), lty = c(1, NA, 1, NA), col = c(col_wv[1], col_ci[1], col_wv[2], col_ci[2]),
cex = 1, pt.cex = c(1.25, 3, 1.25, 3), bty = "n")
}else{
legend(legend_position,
legend = c("Classical WV", "Classical CI", "Robust WV", "Robust CI"),
pch = c(16, 15, 16, 15), lty = c(1, NA, 1, NA), col = c(col_wv[1], col_ci[1], col_wv[2], col_ci[2]),
cex = 1, pt.cex = c(1.25, 3, 1.25, 3), bty = "n")
}
}
}
lines(wv_cl$scales, wv_cl$variance, type = "l", col = col_wv[1], pch = 16)
lines(wv_cl$scales, wv_cl$variance, type = "p", col = col_wv[1], pch = 16, cex = 1.25)
lines(wv_cl$scales, wv_rob$variance, type = "l", col = col_wv[2], pch = 16)
lines(wv_cl$scales, wv_rob$variance, type = "p", col = col_wv[2], pch = 16, cex = 1.25)
}
#' @title Multi-Plot Comparison Between Multiple Wavelet Variances
#' @description
#' This is a helper function for the \code{compare_var()} function.
#' This method accepts the same set of arguments as \code{compare_wvar} and returns a comparision
#' of multiple wavelet variances of different time series accounting for CI values as a set of different plots.
#'
#' @param graph_details List of inputs
#'
#' @author Stephane Guerrier, Justin Lee, and Nathanael Claussen
#' @export
#'
compare_wvar_split = function(graph_details){
old_pars = par(mfrow = c(graph_details$obj_len, graph_details$obj_len),
mar = c(0.5,0.5,0.5,1.5), oma = c(4,4,4,4))
on.exit(par(old_pars))
for (i in 1:graph_details$obj_len){
for (j in 1:graph_details$obj_len){
# Main plot
plot(NA, xlim = graph_details$x_range, ylim = graph_details$y_range, log = "xy", xaxt = 'n',
yaxt = 'n', bty = "n", ann = FALSE)
win_dim = par("usr")
kill_y_tick = graph_details$y_at < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
# Add grid
abline(v = graph_details$x_at, lty = 1, col = "grey95")
abline(h = graph_details$y_at, lty = 1, col = "grey95")
# Add axes and box
box(col = "grey")
# Corner left piece
if (j == 1){
axis(2, at = graph_details$y_at[kill_y_tick],
labels = graph_details$y_labels[kill_y_tick], padj = -0.2, cex.axis = 1/log(graph_details$obj_len))
}
# Corner bottom
if (i == graph_details$obj_len){
axis(1, at = graph_details$x_at, labels = graph_details$x_labels,
padj = 0.1, cex.axis = 1/log(graph_details$obj_len)) # figure out how to size these things for smaller plots
}
# Diag graph
if (i == j){
scales = graph_details$obj_list[[i]]$scales
ci_low = graph_details$obj_list[[i]]$ci_low
ci_high = graph_details$obj_list[[i]]$ci_high
variance = graph_details$obj_list[[i]]$variance
if(graph_details$ci_wv[i] == TRUE){
polygon(c(scales, rev(scales)), c(ci_low, rev(ci_high)),
border = NA, col = graph_details$col_ci[i])
}
lines(scales, variance, type = "l", col = graph_details$col_wv[i], pch = 16)
lines(scales, variance, type = "p", col = graph_details$col_wv[i],
pch = 17, cex = graph_details$point_cex[i]/1.25)
win_dim = par("usr")
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
box()
# if (graph_details$add_legend){
# if (i == j){
# legend(graph_details$legend_position, graph_details$names, bty = "n",
# lwd = 1, pt.cex = graph_details$point_cex, pch = graph_details$point_pch,
# col = graph_details$col_wv, cex=0.7)
# }
# }
}
if (i != j){
scales = graph_details$obj_list[[i]]$scales
ci_low = graph_details$obj_list[[i]]$ci_low
ci_high = graph_details$obj_list[[i]]$ci_high
variance = graph_details$obj_list[[i]]$variance
win_dim = par("usr")
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
if (is.null(graph_details$main[i,j])){
main = paste("WV:", graph_details$names[i],
"vs", graph_details$names[j])
}
box()
if (i < j && graph_details$ci_wv[i] == TRUE){
polygon(c(scales, rev(scales)), c(ci_low, rev(ci_high)),
border = NA, col = graph_details$col_ci[i])
}
lines(scales, variance, type = "l", col = graph_details$col_wv[i], pch = 16)
lines(scales, variance, type = "p", col = graph_details$col_wv[i],
pch = 17, cex = graph_details$point_cex[i]/1.25)
scales = graph_details$obj_list[[j]]$scales
ci_low = graph_details$obj_list[[j]]$ci_low
ci_high = graph_details$obj_list[[j]]$ci_high
variance = graph_details$obj_list[[j]]$variance
if (i < j && graph_details$ci_wv[i] == TRUE){ # don't show confidence intervals
polygon(c(scales, rev(scales)), c(ci_low, rev(ci_high)),
border = NA, col = graph_details$col_ci[j])
}
lines(scales, variance, type = "l", col = graph_details$col_wv[j], pch = 16)
lines(scales, variance, type = "p", col = graph_details$col_wv[j],
pch = 17, cex = graph_details$point_cex[j]/1.25)
}
# Add Details
# @todo: expand win_dim and position $names
if(j==4){
x_vec = 13^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
#mtext(graph_details$names[i], side = 4, line = 0.1, cex = 0.8)
par(xpd = TRUE) #Draw outside plot area
text(x = x_vec[2], y = 0.03, graph_details$names[i], srt = 270, cex = 1.3, col = graph_details$col_wv[i])
par(xpd = FALSE)
}
if(i==1){
mtext(graph_details$names[j], side = 3, line = 0.2, cex = 0.8, col = graph_details$col_wv[j])
}
}
}
mtext(graph_details$ylab, side = 2, line = 2.80, cex = graph_details$cex_labels, outer = T)
mtext(graph_details$xlab, side = 1, line = 2.75, cex = graph_details$cex_labels, outer = T)
if (is.null(graph_details$main)){
main = "Haar Wavelet Variance Representation"
}else{
mtext(main, side = 3, line = 1)
}
}
#' @title Combined Plot Comparison Between Multiple Wavelet Variances
#' @description
#' This is a helper function for the \code{compare_var()} function.
#' This method accepts the same set of arguments as \code{compare_wvar} and returns a single plot
#' that compares multiple wavelet variances of different time series accounting for CI values.
#'
#' @param graph_details List of inputs
#'
#' @author Stephane Guerrier, Justin Lee, and Nathanael Claussen
#' @export
#'
compare_wvar_no_split = function(graph_details){
# Main plot
plot(NA, xlim = graph_details$x_range, ylim = graph_details$y_range, log = "xy", xaxt = 'n',
yaxt = 'n', bty = "n", ann = FALSE)
win_dim = par("usr")
# Main Plot
par(new = TRUE)
plot(NA, xlim = graph_details$x_range, ylim = 10^c(win_dim[3], win_dim[4] + 0.09*(win_dim[4] - win_dim[3])),
log = "xy", xaxt = 'n', yaxt = 'n', bty = "n",
xlab = graph_details$xlab, ylab = graph_details$ylab,
cex.lab = graph_details$cex_labels)
win_dim = par("usr")
# Add Grid
abline(v = graph_details$x_at, lty = 1, col = "grey95")
abline(h = graph_details$y_at, lty = 1, col = "grey95")
# Add Title
x_vec = 10^c(win_dim[1], win_dim[2], win_dim[2], win_dim[1])
y_vec = 10^c(win_dim[4], win_dim[4],
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]),
win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
polygon(x_vec, y_vec, col = "grey95", border = NA)
text(x = 10^mean(c(win_dim[1], win_dim[2])), y = 10^(win_dim[4] - 0.09/2*(win_dim[4] - win_dim[3])),
graph_details$main)
# Add Axes and Box
lines(x_vec[1:2], rep(10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3])),2), col = 1)
y_ticks = graph_details$y_ticks[(10^graph_details$y_ticks) < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))]
kill_y_tick = graph_details$y_at < 10^(win_dim[4] - 0.09*(win_dim[4] - win_dim[3]))
box()
axis(1, at = graph_details$x_at, labels = graph_details$x_labels, padj = 0.3)
axis(2, at = graph_details$y_at[kill_y_tick],
labels = graph_details$y_labels[kill_y_tick], padj = -0.2)
for (i in 1:graph_details$obj_len){
scales = graph_details$obj_list[[i]]$scales
ci_low = graph_details$obj_list[[i]]$ci_low
ci_high = graph_details$obj_list[[i]]$ci_high
variance = graph_details$obj_list[[i]]$variance
if (graph_details$ci_wv[i]){
polygon(c(scales, rev(scales)), c(ci_low, rev(ci_high)),
border = NA, col = graph_details$col_ci[i])
}
lines(scales, variance, type = "l", col = graph_details$col_wv[i], pch = 16)
lines(scales, variance, type = "p", col = graph_details$col_wv[i],
pch = graph_details$point_pch[i], cex = graph_details$point_cex[i])
}
if (graph_details$add_legend){
legend(graph_details$legend_position, graph_details$names, bty = "n",
lwd = 1, pt.cex = graph_details$point_cex, pch = graph_details$point_pch,
col = graph_details$col_wv)
}
}
#' @title Comparison Between Multiple Wavelet Variances
#' @description
#' Displays plots of multiple wavelet variances of different time series accounting for CI values.
#'
#' @param ... One or more time series objects.
#' @param split A \code{boolean} that, if TRUE, arranges the plots into a matrix-like format.
#' @param add_legend A \code{boolean} that, if TRUE, adds a legend to the plot.
#' @param units A \code{string} that specifies the units of time plotted on the x axes. Note: This argument will not be used if xlab is specified.
#' @param xlab A \code{string} that gives a title for the x axes.
#' @param ylab A \code{string} that gives a title for the y axes.
#' @param main A \code{string} that gives an overall title for the plot.
#' @param col_wv A \code{string} that specifies the color of the wavelet variance lines.
#' @param col_ci A \code{string} that specifies the color of the confidence interval shade.
#' @param nb_ticks_x An \code{integer} that specifies the maximum number of ticks for the x-axis.
#' @param nb_ticks_y An \code{integer} that specifies the maximum number of ticks for the y-axis.
#' @param legend_position A \code{string} that specifies the position of the legend (use \code{legend_position = NA} to remove legend).
#' @param ci_wv A \code{boolean} that determines whether confidence interval polygons will be drawn.
#' @param point_cex A \code{double} that specifies the size of each symbol to be plotted.
#' @param point_pch A \code{double} that specifies the symbol type to be plotted.
#' @param names A \code{string} that specifies the name of the WVAR objects.
#' @param cex_labels A \code{double} that specifies the magnification of the labels (x and y).
#' @author Stephane Guerrier and Justin Lee
#' @export
#' @examples
#' set.seed(999)
#' n = 10^4
#' Xt = arima.sim(n = n, list(ar = 0.10))
#' Yt = arima.sim(n = n, list(ar = 0.35))
#' Zt = arima.sim(n = n, list(ar = 0.70))
#' Wt = arima.sim(n = n, list(ar = 0.95))
#'
#' wv_Xt = wvar(Xt)
#' wv_Yt = wvar(Yt)
#' wv_Zt = wvar(Zt)
#' wv_Wt = wvar(Wt)
#'
#' compare_wvar(wv_Xt, wv_Yt, wv_Zt, wv_Wt)
compare_wvar = function(... , split = FALSE, add_legend = TRUE, units = NULL, xlab = NULL,
ylab = NULL, main = NULL, col_wv = NULL, col_ci = NULL, nb_ticks_x = NULL,
nb_ticks_y = NULL, legend_position = NULL, ci_wv = NULL, point_cex = NULL,
point_pch = NULL, names = NULL, cex_labels = 0.8){
obj_list = list(...)
obj_name = as.character(substitute(...()))
obj_len = length(obj_list)
# Check if passed objects are of class wvar
is_wvar = sapply(obj_list, FUN = is, class2 = 'wvar')
if(!all(is_wvar == T)){
stop("Supplied objects must be 'wvar' objects.")
}
# Check length of time series argument
if (obj_len == 0){
stop('No object given!')
}else if (obj_len == 1){
# -> plot.wvar
plot.wvar(..., nb_ticks_x = nb_ticks_x, nb_ticks_y = nb_ticks_y)
}else{
if (is.null(xlab)){
if (is.null(units)){
xlab = expression(paste("Scale ", tau, sep =""))
}else{
xlab = bquote(paste("Scale ", "(", .(units), ")", sep = " "))
}
}else{
xlab = xlab
}
if (is.null(ylab)){
ylab = bquote(paste("Wavelet Variance ", nu^2, sep = " "))
}else{
ylab = ylab
}
if (is.null(ci_wv)){
ci_wv = rep(TRUE, obj_len)
}else{
ci_wv = rep(ci_wv, obj_len)
}
# Main Title
if (split == FALSE){
if (is.null(main)){
main = "Haar Wavelet Variance Representation"
}
}else{
if (!is.null(main) && (dim(main)[1] != obj_len || dim(main)[2] != obj_len)){
main = NULL
}
}
hues = seq(15, 375, length = obj_len + 1)
# Line Colors
if (is.null(col_wv)){
col_wv = hcl(h = hues, l = 65, c = 200, alpha = 1)[seq_len(obj_len)]
}else{
if (length(col_wv) != obj_len){
col_wv = hcl(h = hues, l = 65, c = 200, alpha = 1)[seq_len(obj_len)]
}
}
# CI Colors
if (is.null(col_ci)){
col_ci = hcl(h = hues, l = 80, c = 100, alpha = 0.2)[seq_len(obj_len)]
}else{
if (length(col_ci) != obj_len){
col_ci = hcl(h = hues, l = 80, c = 100, alpha = 0.2)[seq_len(obj_len)]
}
}
# X and Y Limits
x_range = y_range = rep(NULL, 2)
for (i in 1:obj_len){
x_range = range(c(x_range, obj_list[[i]]$scales))
y_range = range(c(y_range, obj_list[[i]]$ci_low, obj_list[[i]]$ci_high))
}
x_low = floor(log10(x_range[1]))
x_high = ceiling(log10(x_range[2]))
y_low = floor(log10(y_range[1]))
y_high = ceiling(log10(y_range[2]))
# Axes Labels and Ticks
if (is.null(nb_ticks_x)){
nb_ticks_x = 6
}
if (is.null(nb_ticks_y)){
nb_ticks_y = 5
}
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > nb_ticks_x){
x_ticks = x_low + ceiling((x_high - x_low)/(nb_ticks_x + 1))*(0:nb_ticks_x)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(10^ .(i))))
x_at = 10^x_ticks
x_actual_length = sum((x_at < x_range[2])*(x_at > x_range[1]))
if (x_actual_length < (3 + as.numeric(split == FALSE))){
x_low = floor(log2(x_range[1]))
x_high = ceiling(log2(x_range[2]))
x_ticks = seq(x_low, x_high, by = 1)
if (length(x_ticks) > 8){
x_ticks = seq(x_low, x_high, by = 2)
}
x_labels = sapply(x_ticks, function(i) as.expression(bquote(2^ .(i))))
x_at = 2^x_ticks
}
y_ticks <- seq(y_low, y_high, by = 1)
if (length(y_ticks) > nb_ticks_y){
y_ticks = y_low + ceiling((y_high - y_low)/(nb_ticks_y + 1))*(0:nb_ticks_y)
}
y_labels = sapply(y_ticks, function(i) as.expression(bquote(10^ .(i))))
y_at = 10^y_ticks
# Legend position
if (is.null(legend_position)){
inter = rep(NA, obj_len)
for (i in 1:obj_len){
inter[i] = obj_list[[i]]$variance[1]
}
mean_wv_1 = mean(inter)
if (which.min(abs(c(y_low, y_high) - log2(mean_wv_1))) == 1){
legend_position = "topleft"
}else{
legend_position = "bottomleft"
}
}
# Type of Points
if (is.null(point_pch)){
inter = rep(15:18, obj_len)
point_pch = inter[1:obj_len]
}else{
if (length(point_pch) != obj_len){
inter = rep(15:18, obj_len)
point_pch = inter[1:obj_len]
}
}
#Size of Points
if (is.null(point_cex)){
inter = rep(c(1.25,1.25,1.25,1.25), obj_len)
point_cex = inter[1:obj_len]
}else{
if (length(point_pch) != obj_len){
inter = rep(c(1.25,1.25,1.25,1.25), obj_len)
point_cex = inter[1:obj_len]
}
}
# Names of WVAR Objects
if (is.null(names)){
names = obj_name
}else{
if (length(names) != obj_len){
names = obj_name
}
}
# Arguments passed into compare_wvar_split or compare_wvar_no_split
graph_details = list(obj_list = obj_list, obj_len = obj_len, names = names, xlab = xlab,
ylab = ylab, col_wv = col_wv, add_legend = add_legend,
col_ci = col_ci, main = main, legend_position = legend_position,
ci_wv = ci_wv, point_cex = point_cex, point_pch = point_pch,
x_range = x_range, y_range = y_range, x_ticks = x_ticks,
x_labels = x_labels, y_labels = y_labels, x_at = x_at, y_at = y_at,
y_ticks = y_ticks, nb_ticks_x = nb_ticks_x, nb_ticks_y = nb_ticks_y,
cex_labels = cex_labels)
if (split == FALSE){
# -> compare_wvar_no_split
compare_wvar_no_split(graph_details)
}else{
# -> compare_wvar_split
compare_wvar_split(graph_details)
}
}
}
Try the wv package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.