Nothing
################################################
## WMTSA package wavelet variance functionality
##
## Functions:
##
## D.statistic
## D.table
## wavVarConfidence
## wavEDOF
##
## Constructor Functions and methods:
##
## wavVar
##
## plot.wavVar
## print.wavVar
## summary.wavVar
## print.summary.wavVar
##
## wavVarTest:
##
## plot.wavVarTest
## print.wavVarTest
##
##################################################
###
# D.statistic
##
"D.statistic" <- function(x)
{
# The D.statistic() function calculates the
# the maximum departure from an expected
# linear increase in cumulative energy
N <- length(x)
Nm1 <- (N - 1)
cumulative.energy <- cumsum(x[1:(length(x)-1)]^2)
P <- cumulative.energy / sum(x^2)
index <- seq(length=Nm1, from=0)
Dplus <- max((index + 1) / Nm1 - P)
Dminus <- max(P - (index / Nm1))
D <- max(c(Dplus, Dminus))
# my (dbp) very first comment!!! This was a great comment!
D
}
###
# D.table
##
"D.table" <- function(n.sample=c(127, 130), significance=c(0.1,0.05,0.01),
lookup=TRUE, n.realization=10000, n.repetition=3, tolerance=1e-6)
{
# define local functions
"D.table.preprocess" <- function(significance, n.sample)
{
# The D.table.preprocess() function sorts the significance
# and n.sample vectors and checks for appropriate ranges.
# The new n.sample and significance vectors are returned
# in a list object.
# test to make sure the significance is in the proper range
significance.range <- range(significance)
checkRange(significance.range)
n.sample <- as.integer(n.sample)
if (!(all(n.sample > 0)))
stop("n.sample value(s) must be positive integers")
# sort the significances from high to low
significance <- unique(rev(sort(significance)))
# sort the n.sample vector from low to high
n.sample <- unique(sort(n.sample))
list(significance=significance, n.sample=n.sample)
}
assign("D.table.preprocess", D.table.preprocess)
"D.table.names" <- function(n.sample, significance)
{
# The D.table.names() function creates the names for the
# rows and columns of a D.table matrix
ilow <- which(n.sample < 128)
method <- rep("Inclan-Tiao",length(n.sample))
method[ilow] <- rep("Monte Carlo",length(ilow))
list(method , c("N", paste(significance*100,"%",sep="")))
}
assign("D.table.names", D.table.names)
"significance.percentage.to.number" <- function(x)
{
# The significance.percentage.to.number() function
# converts the significance(s) expressed as a
# character vector with a percentage to a number(s)
# For example, "10%" goes to 10.
# get rid of "%" characters and turn column
# names into significances
as.numeric(gsub("%","",x))/100
}
"D.table.make" <- function(significance=c(0.1, 0.05, 0.01), n.sample=c(10,200),
n.repetition=3, n.realization=10000, tolerance=1e-6)
{
# The D.table.make() function returns a matrix of D-statistics
# corresponding to a single round of Monte Carlo simulations
# if n.sample < 128 combined with the D-statistics
# calculated via Incan-Tiao for n.sample >= 128.
# The rows are the realization lengths while the columns are
# the significances.
# define local functions
"D.MonteCarlo" <- function(n.sample=10, n.realization=10000, significance=0.1)
{
# The D.MonteCarlo() function returns the D-statistics corresponding to
# the significance(s) for a single realization length (n.sample < 128)
# precalculate the indices of the sorted D-statistic vector
# corresponding to the significances
index <- floor((1 - significance) * n.realization)
D.realization <- sort(unlist(lapply(seq(length=n.realization),
function(i, n.sample){
# generate a N(0,1) realization
G <- rnorm(n=n.sample, mean=0, sd=1)
# calculate and return the D-statistic
D.statistic(G)
},
n.sample=n.sample)))
# return D-statistics corresponding to the significances
(D.realization[index] + D.realization[index + 1]) / 2
}
"D.InclanTiao" <- function(n.sample=10, tolerance=1e-6, significance=0.1)
{
# define local functions
"InclanTiao.root" <- function(d, n.sample=100, tolerance=1e-6, significance=0.1)
{
# The Inclan-Tiao root function is used in
# combination with the uniroot function to approximate
# the D-statistic for sample sizes larger than (say) 128.
# The true function involves an infinite summation of exponentials
# which tend to decay rather quickly. Thus, a tolerance
# denoting the amplitude threshold limits the number of
# of summations and increases the speed of calculation.
if (n.sample > 100){
max.iterations <- ceiling(sqrt(-log(tolerance/2) / n.sample / d^2))
l <- seq(length=max.iterations)
2 * sum((-1)^(l+1) * exp(-l*l*n.sample*d*d)) - significance
}
else
stop("The Inclan-Tiao method works best for N.SAMPLE > 100")
}
# returns the D-statistics corresponding to
# the significance(s) for a single realization length (n.sample >= 128)
unlist(lapply(seq(along=significance),
function(i, n.sample, significance, tolerance,InclanTiao.root){
eps <- .Machine$double.eps
uniroot(InclanTiao.root,
c(1e-4,1-eps), n.sample=n.sample,
significance=significance[i],
tolerance=tolerance)$root
},
n.sample=n.sample,
significance=significance,
tolerance=tolerance,
InclanTiao.root=InclanTiao.root))
}
# PRIMARY CODE SET
# check ranges on input and sort
input <- D.table.preprocess(significance, n.sample)
significance <- input$significance
n.sample <- input$n.sample
D.nsample <- unlist(lapply(seq(along=n.sample),
function(i, n.sample, n.realization, significance, tolerance,
n.repetition, D.MonteCarlo, D.InclanTiao){
# calculate the critical levels corresponding to
# each realization length (n.sample)
if (n.sample[i] < 128){
D.round <- lapply(seq(length=n.repetition),
function(j, n.sample, n.realization, significance, D.MonteCarlo){
D.MonteCarlo(n.sample=n.sample,
n.realization=n.realization,
significance=significance)
},
n.sample=n.sample[i],
n.realization=n.realization,
significance=significance,
D.MonteCarlo=D.MonteCarlo)
D.array <- array(unlist(D.round), c(length(significance), n.repetition))
D <- rowMeans(D.array)
}
else{
D <- D.InclanTiao(n.sample=n.sample[i],
tolerance=tolerance,
significance=significance)
}
return(D)
},
n.sample=n.sample,
n.realization=n.realization,
significance=significance,
tolerance=tolerance,
n.repetition=n.repetition,
D.MonteCarlo=D.MonteCarlo,
D.InclanTiao=D.InclanTiao))
D.tab <- cbind(n.sample, matrix(D.nsample, ncol=length(significance), byrow=TRUE))
dimnames(D.tab) <- D.table.names(n.sample, significance)
return(D.tab)
}
"D.table.fill.missing" <- function(D.row, D.table.make, col.names=NULL)
{
# The D.table.fill.missing() function fills in the
# the D-statistics missing from a D.table row
# define local functions
"significance.percentage.to.number" <- function(x)
{
# The significance.percentage.to.number() function
# converts the significance(s) expressed as a
# character vector with a percentage to a number(s)
# For example, "10%" goes to 10.
# get rid of "%" characters
x <- gsub("%","",x)
# ... turn them into significances
as.numeric(x)/100
}
missing.indices <- is.na(D.row)
if (any(missing.indices)){
significance.na <- significance.percentage.to.number(col.names[missing.indices])
n.sample.na <- D.row[1]
D.na <- D.table.make(significance=significance.na,
n.sample=n.sample.na,
n.realization=100)
# invoke column offset for significances returned by the
# D.table.make() function (the first column denotes the n.sample
# value, so we skip this column when combining results)
D.na.indices <- seq(from=2,length=length(which(missing.indices)))
D.row[missing.indices] <- D.na[D.na.indices]
}
D.row
}
if (lookup){
# check ranges on input and sort
input <- D.table.preprocess(significance, n.sample)
significance <- input$significance
n.sample <- input$n.sample
# verify existence of critical D statistic table.
# if not available, issue error
if (!exists("D.table.critical")){
warning("The D.table.critical lookup table does not exist. The table is calculated explicitly")
D <- D.table.make(significance=significance,
n.sample=n.sample,
n.repetition=n.repetition,
n.realization=n.realization,
tolerance=tolerance)
}
else{
col.names <- dimnames(D.table.critical)[[2]]
# find column names that have a "%" in them
ipercentage <- grep("%", col.names)
col.names <- col.names[ipercentage]
# ... turn them into significances
table.significance <- significance.percentage.to.number(col.names)
# find the column indices in D.critical which match the
# significances. Non-matching elements of the significance
# vector are returned as NA's
column.offset <- min(ipercentage) - 1
cols <- match(significance, table.significance) + column.offset
# find the row indices in in D.critical which match the
# n.sample values. Non-matching elements of the n.sample
# vector are returned as NA's
rows <- match(n.sample, D.table.critical[,1])
# build critical D-statistic table
D <- cbind(N=n.sample, D.table.critical[rows,cols, drop=FALSE])
dimnames(D) <- D.table.names(n.sample, significance)
# if there are missing values in the table, fill them in
if (any(is.na(D))){
D <- t(apply(D, MARGIN=1,
function(x, D.table.fill.missing, D.table.make, col.names){
D.table.fill.missing(x, D.table.make, col.names)
},
D.table.fill.missing=D.table.fill.missing,
D.table.make=D.table.make,
col.names=dimnames(D)[[2]]))
}
} # end check for existence of D.table.critical
} # end if lookup
else{
D <- D.table.make(significance=significance,
n.sample=n.sample,
n.repetition=n.repetition,
n.realization=n.realization,
tolerance=tolerance)
}
dimnames(D) <- D.table.names(n.sample, significance)
# remove local functions
remove("D.table.preprocess")
remove("D.table.names")
return(D)
}
###
# wavVarConfidence
###
"wavVarConfidence" <- function(wvar, edof, probability=0.95)
{
# turn EDOF list into a vector, one EDOF element per scale
edof <- unlist(edof)
if (any(is.na(edof))){
na.conf <- rep(NA, length(edof))
confidence <- list(low=na.conf, high=na.conf)
}
else{
confidence <- itCall("RS_wavelets_variance_confidence",
wvar,
edof,
probability)
#
#COPY=rep(FALSE,3),
#CLASSES=c("numeric","numeric","numeric"),
#PACKAGE="ifultools")
confidence <- lapply(confidence,function(x,crystals){
y <- as.vector(x)
names(y) <- crystals
return(y)
},
crystals=names(edof))
names(confidence) <- c("low","high")
}
return(confidence)
}
###
# wavEDOF
###
"wavEDOF" <- function(x, wavelet="s8", levels=NULL, sdf=NULL, sdfargs=NULL, sampling.interval=1, n.fft=1024)
{
if (is(x,"wavTransform")){
# obtain interior wavelet coefficients
Wb <- wavBoundary(x)
J <- x$dictionary$n.levels
# form levels vector if empty
if (is.null(levels))
levels <- seq(J)
if (max(levels) > J)
stop("Levels exceeds maximum in transform")
crystals <- paste("d", levels, sep="")
interior <- Wb$interior[crystals]
n.coeff <- Wb$interior$length[crystals]
# form the unbiased wavelet variance estimates
variance <- unlist(lapply(interior,function(x) sum(x^2))) / n.coeff
# adjust normalization for DWT style
if (x$xform == "dwt")
variance <- variance / 2^(levels)
# map filters
filter <- mutilsFilterType(x$dictionary$wavelet)
# evaluate input SDF over frequencies [0, Nyquist]
Sx <- mutilsSDF(sdf, sdfargs=sdfargs, n.freq=n.fft, sampling.interval=sampling.interval)
# adjust for case where there are no interior coefficients
# at a given level
bad <- which(is.na(variance))
if (length(bad)){
interior <- interior[-bad]
n.coeff <- n.coeff[-bad]
levels <- levels[levels != bad]
crystals <- crystals[-bad]
variance <- variance[-bad]
}
# calculate the EDOF
edof <- itCall("RS_wavelets_variance_edof",
as.numeric(unlist(interior)),
as.integer(n.coeff),
as.numeric(variance),
as.integer(levels),
as.numeric(Sx),
as.integer(filter$type),
as.integer(filter$length))
#
#COPY=rep(FALSE,7),
#CLASSES=c("numeric","integer","numeric","integer","numeric","integer","integer"),
#PACKAGE="ifultools")
edof <- lapply(edof,
function(x,crystals){
y <- as.vector(x)
names(y) <- crystals
invisible(y)}, crystals=crystals)
names(edof) <- paste("EDOF", 1:3, sep="")
names(n.coeff) <- crystals
# pack results into a list
edof <- c(edof, list(variance.unbiased=variance, n.coeff=n.coeff))
# if any members of edof are < 0, that means that they
# were not calculated. in this case, place an NA for
# each element of that edof mode
if (any (edof[[2]] < 0))
edof[[2]] <- rep(NA,length(edof[[2]]))
return(edof)
}
else if (is.element(class(x), c("numeric","ts","rts","signalSeries"))){
y <- create.signalSeries(x)
dt <- y@positions@by
if (is.null(levels))
n.level <- ilogb(length(y),base=2)
else
n.level <- as.integer(max(levels))
xform <- wavMODWT(y, wavelet=wavelet, n.levels=n.level)
return(wavEDOF(xform, sdf=sdf, levels=levels, sampling.interval=dt, sdfargs=sdfargs))
}
else{
stop("Class of input is unsupported")
}
}
##############################################################################
##
## Constructor wavVar:
##
## plot.wavVar
## print.wavVar
## summary.wavVar
## print.summary.wavVar
##
##############################################################################
###
# Constructor: wavVar
###
"wavVar" <- function(x, xform="modwt", wavelet="s8", n.levels=NULL,
position=list(from=1,by=1,units=character()), units=character(),
documentation=character(), sdf=NULL, sdfargs=NULL, sampling.interval=1, n.fft=1024)
{
# get the series name from the parent
# (or from this function if this is called as a top level function)
series.name <- deparse(substitute(x))
series <- create.signalSeries(x, position=position, units=units,
title.data=series.name, documentation=documentation)
# map filters
filter <- mutilsFilterType(wavelet)
if (is.null(n.levels))
n.levels <- wavMaxLevel(n.taps=filter$length, n.sample=length(x), xform=xform)
# evaluate input SDF over frequencies [0, Nyquist]
Sx <- mutilsSDF(sdf=sdf, sdfargs=sdfargs, n.freq=n.fft, sampling.interval=deltat(x))
type <- mutilsTransformType(xform)
obj <- itCall("RS_wavelets_variance",
as.numeric(series@data),
as.integer(type),
as.integer(filter$type),
as.integer(filter$length),
as.integer(n.levels),
as.numeric(Sx))
#
#COPY=rep(FALSE,6),
#CLASSES=c("numeric","integer","integer","integer","integer","numeric"),
#PACKAGE="ifultools")
is.conf <- as.logical(lowerCase(xform) == "modwt")
# develop crystal names
levels <- seq(along=obj[[2]][[1]])
crystals <- paste("d",levels,sep="")
sampling.interval <- series@positions@by
scales <- 2^(levels - 1)*sampling.interval
# name the data
names(scales) <- crystals
vtime <- lapply(obj[[1]], as.vector)
names(vtime) <- c("biased","unbiased")
vblock <- lapply(obj[[2]], function(x,crystals){
y <- as.vector(x);names(y) <- crystals;return(y) }, crystals=crystals)
names(vblock) <- c("biased","unbiased")
filters <- wavDaubechies(wavelet=wavelet, normalized=(xform == "modwt"))
if (is.conf){
conf <- lapply(obj[[3]], function(x,crystals){
low <- x[1,]
high <- x[2,]
names(low) <- crystals
names(high) <- crystals
return(list(low=low,high=high))
},
crystals=crystals)
names(conf) <- paste("n",1:3,sep="")
separation <- itCall("RS_wavelets_transform_coefficient_boundaries",
as.integer(n.levels), as.integer(filters$length), as.integer(length(series)), as.integer(type))
#
#COPY=rep(FALSE,4),
#CLASSES=rep("integer",4),
#PACKAGE="ifultools")
unbiased.length <- separation[[3]][levels]
names(unbiased.length) <- crystals
edof <- lapply(obj[[4]],
function(x,crystals){
y <- as.vector(x)
names(y) <- crystals
return(y)
}, crystals=crystals)
names(edof) <- paste("EDOF",1:3,sep="")
}
else{
conf <- list(EDOF1=NA, EDOF2=NA, EDOF3=NA)
edof <- list("n1"=NA, "n2"=NA, "n3"=NA)
unbiased.length <- NA
}
# create dictionary
dict <- list(wavelet =wavelet,
dual = FALSE,
decimate = FALSE,
n.sample = length(x),
attr.x = NULL,
n.levels = n.levels,
boundary = "periodic",
conv = TRUE,
analysis.filter=list(low=filters$scaling, high=filters$wavelet),
synthesis.filter=list(low=filters$scaling, high=filters$wavelet),
fast=TRUE,
is.complex=FALSE)
oldClass(dict) <- "wavDictionary"
series@title <- series.name
wvar <- list(time=vtime,
block=vblock,
scales=scales,
xform=xform,
dictionary=dict,
series=series,
sampling.interval=sampling.interval)
z <- c(wvar, list(confidence=c(edof,conf,list(length=unbiased.length))))
if (is.null(sdf))
z$confidence$EDOF2 <- NA
oldClass(z) <- "wavVar"
z
}
###
# plot.wavVar
###
"plot.wavVar" <- function(x, type="unbiased", time=FALSE, xlab=NULL,
title=NULL, ylab=NULL, edof=3, pch=18, logfunc=log10, ...)
{
# check input data
checkVectorType(edof,"integer")
checkRange(edof,c(1,3))
# now that the internal functions are developed for plotting,
# begin plot method code ...
# map the plot type
biased <- type == "biased"
# obtain sequence name and, if too long, truncate
series <- x$series
series.name <- series@title
sampling.interval <- x$sampling.interval
if (length(series.name) < 1)
series.name <- "x"
if (length(sampling.interval) < 1)
sampling.interval <- 1
# obtain transform type, and wavelet type
xform <- upperCase(x$xform)
wavelet <- x$dictionary$wavelet
wavelet.str <- paste("(", wavelet, ")")
biased.title <- paste("Biased", xform, "Wavelet Variance", wavelet.str)
unbiased.title <- paste("Unbiased", xform, "Wavelet Variance", wavelet.str)
# create matrix of variance data
biased.data <- x$block$biased
unbiased.data <- x$block$unbiased
# obtain the crystals for each type of data
unbiased.crystals <- names(unbiased.data)
biased.crystals <- names(biased.data)
if (is.null(xlab)){
if (time == FALSE) {
xlab <- "Scale"
xunits <- x$series@units.position
if (length(xunits) > 0)
xlab <- paste(xlab, " (", xunits, ")",sep="")
}
else
xlab <- "Time"
}
if (is.null(title))
title <- paste("Wavelet Variance of", series.name)
if (time){
# time-dependent wavelet variance plots
old.plt <- splitplot(1,1,1)
on.exit(par(old.plt))
same.scale <- TRUE # make sure magnitudes are similar
zero.line <- TRUE
N <- x$dictionary$n.sample
dt <- x$series@positions@by
times <- (1:N) * dt
nms <- rev(names(x$time$unbiased))
if (biased){
wavStackPlot(x$time$biased[nms], times=times, same.scale=same.scale, zeroline=zero.line)
title(paste(biased.title,"\n\n",series.name),xlab=xlab)
}
else{
wavStackPlot(x$time$unbiased[nms], times=times, same.scale=same.scale, zeroline=zero.line)
title(paste(unbiased.title,"\n\n",series.name),xlab=xlab)
}
}
else{
# block-dependent wavelet variance plot
xdata <- logfunc(x$scales)
ydata <- logfunc(ifelse1(biased, biased.data, unbiased.data))
if (any(is.na(x$confidence$EDOF2)))
edof <- edof[edof != 2]
if (!length(edof))
stop("No data exists for current EDOF specification")
conf <- ifelse1(xform == "DWT", list(), x$confidence[paste("n", edof, sep="")])
# use the same plot limits for all EDOF plots for comparison purposes
xlim <- range(xdata, na.rm=TRUE)
if (length(conf))
ylim <- range(logfunc(unlist(conf)), na.rm=TRUE)
logtext <- deparse(substitute(logfunc))
xlab <- paste(logtext, "(scale)", sep="")
ylab <- paste(logtext, ifelse1(biased, biased.title, unbiased.title))
for (i in seq(length(edof))){
old.plt <- splitplot(1,length(edof),i, gap=0.1)
if (i == 1)
on.exit(par(old.plt))
if (length(conf)){
low <- logfunc(conf[[i]]$low)
high <- logfunc(conf[[i]]$high)
plot(xlim, ylim, type="n", xlab=xlab, ylab=ifelse1(i==1,ylab,""))
arrows(xdata, low, xdata, high, code=3, col="black", angle=90, lwd=2, length=0.1)
points(xdata, ydata, pch=19, col="blue", cex=1.5)
mtext(paste("EDOF", edof[i]), side=3, adj=1, line=0.5)
}
else
plot(xdata, ydata, pch=19, col="blue", cex=1.5)
grid()
}
}
invisible(NULL)
}
###
# print.wavVar
###
"print.wavVar" <- function(x, justify="left", sep=":", ...)
{
series <- x$series
title <- series@title
sampling.interval <- x$sampling.interval
pos.units <- series@units.position
if (length(title) < 1)
title <- "x"
crystals <- wavSortCrystals(names(x$block$biased))
dict <- x$dict
main <- paste(upperCase(x$xform)," wavelet variance of ", title, sep="")
z <- c(as.list(x$dict), list("Sampling interval"=sampling.interval,
"Units"=ifelse1(length(pos.units) > 1, paste("(", pos.units, ")"), NULL)))
prettyPrintList(as.list(x$dict), header=main, sep=sep, justify=justify, ...)
scales <- x$scales[crystals]
print(matrix(scales,nrow=1,dimnames=list("Scales",crystals)), digits=2)
invisible(x)
}
###
# summary.wavVar
###
"summary.wavVar" <- function(object, ...)
{
# create matrix of variance data
biased.data <- object$block$biased
unbiased.data <- object$block$unbiased
# obtain the crystals for each type of data
unbiased.crystals <- wavSortCrystals(names(unbiased.data))
biased.crystals <- wavSortCrystals(names(biased.data ))
# bind the variance data into a row major matrix.
# Since the unbiased estimates use only interior
# wavelet coefficients, and since some crystals
# may not contain any interior coefficients,
# those crystals are ignored. Consequently, the
# length of the biased and unbiased wavelet vectors
# may differ.
data <- matrix(NA,nrow=2,ncol=length(biased.data))
dimnames(data) <- list(c("Biased","Unbiased"),biased.crystals)
data[1,] <- biased.data[biased.crystals]
data[2,][unbiased.crystals] <- unbiased.data[unbiased.crystals]
# obtain sequence name and, if too long, truncate
series <- object$series
series.name <- series@title
sampling.interval <- object$sampling.interval
if (length(series.name) < 1)
series.name <- "x"
if (length(sampling.interval) < 1)
sampling.interval <- 1
cut <- 15
if (nchar(series.name) > cut)
series.name <- paste(substring(series.name,1,cut), " ... ",sep="")
J <- object$dictionary$n.levels
if (object$xform == "modwt"){
confidence <- object$confidence
# produce matrix to hold data
scale.data <- unlist(object$scales[unbiased.crystals])
EDOF1.data <- unlist(round(confidence$EDOF1[unbiased.crystals]))
EDOF2.data <- unlist(round(confidence$EDOF2[unbiased.crystals]))
EDOF3.data <- unlist(round(confidence$EDOF3[unbiased.crystals]))
length.data <- unlist(confidence$length[unbiased.crystals])
level.data <- ilogb(scale.data/sampling.interval, base=2) + 1
conf.mat <- rbind(level.data, scale.data,
EDOF1.data, EDOF2.data,EDOF3.data,
length.data)
# supply names for each dimension of matrix
dimnames(conf.mat) <- list(c("Level","Scale","[EDOF1]","[EDOF2]","[EDOF3]","# coeffs"),
biased.crystals[1:ncol(conf.mat)])
}
else
conf.mat <- NA
obj <- list(conf.mat=conf.mat,
xform=object$xform,
vmat=data,
series.name=series.name)
oldClass(obj) <- "summary.wavVar"
return(obj)
}
###
# print.summary.wavVar
###
"print.summary.wavVar" <- function(x, digits=max(2, .Options$digits - 4), ...)
{
main <- paste(upperCase(x$xform)," wavelet variance of ", x$series.name,sep="")
cat(main,"\n", rep("-",nchar(main)),"\n",sep="")
print(x$vmat, digits=digits)
if (x$xform == "modwt"){
cat("\nConfidence Interval Data for Unbiased\nTime Independent Wavelet Variance:\n\n")
print(x$conf.mat)
}
invisible(x)
}
##############################################################################
##
## Constructor wavVarTest:
##
## plot.wavVarTest
## print.wavVarTest
##
##############################################################################
###
# wavVarTest (constructor)
##
"wavVarTest" <- function(x, wavelet="s8", n.levels=NULL,
significance=c(0.1,0.05,0.01), lookup=TRUE, n.realization=10000,
n.repetition=3, tolerance=1e-6)
{
# define local functions
"is.dyadic" <- function(x) as.logical(2^ilogb(x,base=2) == x && x > 0)
"D.wavebound" <- function(x, good.crystals)
{
# Calculate the D statistic, which estimates the largest departure from an assumed
# linear increase in cumulative (interior wavelet coefficient)
# energy. If the wavelet energy is homogeneously distributed across time,
# its cumulative energy should increase approximately linearly.
if (!is(x,"wavBoundary"))
stop("Input must be of class \"wavBoundary\"")
# obtain energy interior coefficients
# in each crystal
energy <- x$interior$energy[good.crystals]
# obtain nonzero energy components and corresponding
# crystal names
energy <- energy[energy > 0]
crystals <- names(energy)
# return D-statistic
unlist(lapply(x$interior[crystals], D.statistic))
}
if (is(x,"wavBoundary")){
# check for convolution style filtering
if (!x$dictionary$conv)
stop("Convolution style filtering must be used in the DWT ",
"for for homogeneity of variance tests")
series.name <- wavTitle(x)
}
else{
# get the series name
series.name <- deparseText(substitute(x))
# check for power of 2
dyad <- is.dyadic(length(x))
# perform a dwt on the data
n.sample <- length(x)
if (is.null(n.levels))
n.levels <- wavMaxLevel(
n.taps= mutilsFilterType(wavelet)$length,
n.sample=n.sample,
xform="dwt")
# perform the dwt using convolution style filtering
W <- wavDWT(x, wavelet=wavelet, n.levels=n.levels)
# perform a wavelet boundary separation
WW <- wavBoundary(W)
}
# obtain dwt specifications
dict <- WW$dictionary
n.levels <- dict$n.levels
series <- WW$series
wave.class <- class(W)
crystals <- WW$crystals
# estimate departure from an assumed linear increase in cumulative
# (interior wavelet coefficient) energy
Nj <- WW$interior$length
Nj <- Nj[Nj > 0]
good.crystals <- names(Nj)
# get rid of scaling coefficient crystals
bad <- charmatch("s", good.crystals)
if (!is.na(bad)){
Nj <- Nj[-bad]
good.crystals <- good.crystals[-bad]
}
Dj <- D.wavebound(WW, good.crystals)
# create a D-statistics table
D.critical <- D.table(n.sample=Nj,
significance=significance,
lookup=lookup,
tolerance=tolerance,
n.repetition=n.repetition,
n.realization=n.realization)
# if duplicate Nj exist, replicate that row in D.critical.
# this will occur with sJ and dJ wavelet crystals for example
D.critical <- D.critical[pmatch(Nj,D.critical[,1],duplicates.ok=TRUE),]
# perform the homogeneity test
is.homogeneous <- apply(D.critical[,-1], 2,
function(x,Dj) Dj < x ,Dj=Dj[good.crystals])
dimnames(is.homogeneous)[[2]] <-
dimnames(D.critical)[[2]][-1]
crystal <- names(Dj)
# explicitly set the series name
series@title <- series.name
z <- list(pass=is.homogeneous,
D=Dj,
D.critical=D.critical,
significance=significance,
n.sample=Nj,
dictionary=dict,
series=series,
lookup=lookup,
tolerance=tolerance,
n.repetition=n.repetition,
n.realization=n.realization)
oldClass(z) <- "wavVarTest"
return(z)
}
###
# plot.wavVarTest
##
"plot.wavVarTest" <- function(x, ...)
{
series.name <- wavTitle(x)
sigcols <- 2:ncol(x$D.critical)
key.names <- dimnames(x$D.critical)[[2]][sigcols]
xrange <- range(x$D.critical[,1])
yrange <- range(as.vector(x$D.critical[,2:ncol(x$D.critical)]))
for (i in 1:sigcols[length(sigcols)]){
if (i == 1)
plot(x$n.sample, x$D.critical[,i], ylab="D-statistic",
xlab="Number of Samples", xlim=xrange,ylim=yrange,log="x")
else
points(x$n.sample, x$D.critical[,i], pch=i)
}
# obtain axis limits [x1 x2 y1 y2]
ax <- par("usr")
# add a title
title(series.name)
}
###
# print.wavVarTest
##
"print.wavVarTest" <- function(x, ...)
{
series.name <- wavTitle(x)
dict <- x$dictionary
cat("Homogeneity Test for Discrete Wavelet Transform of", series.name, "\n")
cat("\nPass:\n")
print(x$pass)
cat("\nD-statistic critical values comparison:\n")
methods <- dimnames(x$D.critical)[[1]]
D.cat <- cbind(N=x$D.critical[,1], "D"=x$D, x$D.critical[,-1])
dimnames(D.cat)[[1]] <- paste(dimnames(x$pass)[[1]], methods, sep=": ")
print(D.cat)
cat("\nLookup:", x$lookup,"\n")
if (any(dimnames(x$D.critical)[[1]] == "Inclan-Tiao"))
cat("Inclan-Tiao tolerance: ", x$tolerance,"\n")
if (x$lookup){
cat("Repetitions: 3 (lookup),", x$n.repetition, "(non-lookup)\n")
cat("Realizations/repetition: 10000 (lookup),", x$n.realization, "(non-lookup)\n")
}
else{
cat("Repetitions:", x$n.repetition, "\n")
cat("Realizations/repetition:", x$n.realization, "\n")
}
cat("\n")
print(dict)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.