R: Threshold (DWT) wavelet decomposition object

threshold.wd {wavethresh}

R Documentation

Threshold (DWT) wavelet decomposition object

Description

This function provides various ways to threshold a wd class object.

Usage

## S3 method for class 'wd'
threshold(wd, levels = 3:(nlevelsWT(wd) - 1), type = "soft", policy = "sure",
        by.level = FALSE, value = 0, dev = madmad, boundary = FALSE, verbose = FALSE,
        return.threshold = FALSE, force.sure = FALSE, cvtol = 0.01,
	cvmaxits=500, Q = 0.05, OP1alpha = 0.05, 
        alpha = 0.5, beta = 1, C1 = NA, C2 = NA, C1.start = 100,
	al.check=TRUE, ...)

Arguments

`wd`	The DWT wavelet decomposition object that you wish to threshold.
`levels`	a vector of integers which determines which scale levels are thresholded in the decomposition. Each integer in the vector must refer to a valid level in the `wd` object supplied. This is usually any integer from 0 to `nlevelsWT`(wd)-1 inclusive. Only the levels in this vector contribute to the computation of the threshold and its application.
`type`	determines the type of thresholding this can be "hard" or "soft".
`policy`	selects the technique by which the threshold value is selected. Each policy corresponds to a method in the literature. At present the different policies are: "`universal`", "`LSuniversal`", "`sure`", "`BayesThresh`", "`cv`", "`fdr`", "`op1`", "`op2`", "`manual`", "`mannum`" and "`probability`". The policies are described in detail below.
`by.level`	If FALSE then a global threshold is computed on and applied to all scale levels defined in `levels`. If TRUE a threshold is computed and applied separately to each scale level.
`value`	This argument conveys the user supplied threshold. If the `policy="manual"` then `value` is the actual threshold value; if the `policy="mannum"` then `value` conveys the total number of ordered coefficients kept (from the largest); if `policy="probability"` then `value` conveys the the user supplied quantile level.
`dev`	this argument supplies the function to be used to compute the spread of the absolute values coefficients. The function supplied must return a value of spread on the variance scale (i.e. not standard deviation) such as the `var()` function. A popular, useful and robust alternative is the `madmad` function.
`boundary`	If this argument is TRUE then the boundary bookeeping values are included for thresholding, otherwise they are not.
`verbose`	if TRUE then the function prints out informative messages as it progresses.
`return.threshold`	If this option is TRUE then the actual value of the threshold is returned. If this option is FALSE then a thresholded version of the input is returned.
`force.sure`	If TRUE then the `sure` threshold is computed on a vector even when that vector is very sparse. If FALSE then the normal SUREshrink procedure is followed whereby the universal threshold is used for sparse vectors of coefficients.
`cvtol`	Parameter for the cross-validation `"cv"` policy.
`cvmaxits`	Maximum number of iterations allowed for the cross-validation `"cv"` policy.
`Q`	Parameter for the false discovery rate `"fdr"` policy.
`OP1alpha`	Parameter for Ogden and Parzen's first "`op1`" and `"op2"` policies.
`alpha`	Parameter for BayesThresh `"BayesThresh"` policy.
`beta`	Parameter for BayesThresh `"BayesThresh"` policy.
`C1`	Parameter for BayesThresh `"BayesThresh"` policy.
`C2`	Parameter for BayesThresh `"BayesThresh"` policy.
`C1.start`	Parameter for BayesThresh `"BayesThresh"` policy.
`al.check`	If TRUE then the function checks that the levels are in ascending order. If they are not then this can be an indication that the default level arguments are not appropriate for this data set (`wd` object). However, a strange order might be appropriate for some reason if deliberately set, so setting this argument equal to FALSE turns off the check and warning.
`...`	any other arguments

Details

This function thresholds or shrinks wavelet coefficients stored in a wd object and returns the coefficients in a modified wd object. See the seminal papers by Donoho and Johnstone for explanations about thresholding. For a gentle introduction to wavelet thresholding (or shrinkage as it is sometimes called) see Nason and Silverman, 1994. For more details on each technique see the descriptions of each method below

The basic idea of thresholding is very simple. In a signal plus noise model the wavelet transform of signal is very sparse, the wavelet transform of noise is not (in particular, if the noise is iid Gaussian then so if the noise contained in the wavelet coefficients). Thus since the signal gets concentrated in the wavelet coefficients and the noise remains "spread" out it is "easy" to separate the signal from noise by keeping large coefficients (which correspond to signal) and delete the small ones (which correspond to noise). However, one has to have some idea of the noise level (computed using the dev option in threshold functions). If the noise level is very large then it is possible, as usual, that no signal "sticks up" above the noise.

There are many components to a successful thresholding procedure. Some components have a larger effect than others but the effect is not the same in all practical data situations. Here we give some rough practical guidance, although you must refer to the papers below when using a particular technique. You cannot expect to get excellent performance on all signals unless you fully understand the rationale and limitations of each method below. I am not in favour of the "black-box" approach. The thresholding functions of WaveThresh3 are not a black box: experience and judgement are required!

Some issues to watch for:

levels

The default of levels = 3:(wd$nlevelsWT - 1) for the levels option most certainly does not work globally for all data problems and situations. The level at which thresholding begins (i.e. the given threshold and finer scale wavelets) is called the primary resolution and is unique to a particular problem. In some ways choice of the primary resolution is very similar to choosing the bandwidth in kernel regression albeit on a logarithmic scale. See Hall and Patil, (1995) and Hall and Nason (1997) for more information. For each data problem you need to work out which is the best primary resolution. This can be done by gaining experience at what works best, or using prior knowledge. It is possible to "automatically" choose a "best" primary resolution using cross-validation (but not in WaveThresh).

Secondly the levels argument computes and applies the threshold at the levels specified in the levels argument. It does this for all the levels specified. Sometimes, in wavelet shrinkage, the threshold is computed using only the finest scale coefficients (or more precisely the estimate of the overall noise level). If you want your threshold variance estimate only to use the finest scale coefficients (e.g. with universal thresholding) then you will have to apply the threshold.wd function twice. Once (with levels set equal to nlevelsWT(wd)-1 and with return.threshold=TRUE to return the threshold computed on the finest scale and then apply the threshold function with the manual option supplying the value of the previously computed threshold as the value options.

Thirdly, if you apply wavelet shrinkage to a small data set then you need to ensure you've chosen the levels argument appropriately. For example, if your original data was of length 8, then the associated wd wavelet decomposition object will only have levels 0, 1 and 2. So, the default argument for levels (starting at 3 and higher) will almost certainly be wrong. The code now warns for these situations.

by.level

for a wd object which has come from data with noise that is correlated then you should have a threshold computed for each resolution level. See the paper by Johnstone and Silverman, 1997.

Value

An object of class wd. This object contains the thresholded wavelet coefficients. Note that if the return.threshold option is set to TRUE then the threshold values will be returned rather than the thresholded object.

RELEASE

Note

POLICIES This section gives a brief description of the different thresholding policies available. For further details see the associated papers. If there is no paper available then a small description is provided here. More than one policy may be good for problem, so experiment! They are arranged here in alphabetical order:

BayesThresh: See Abramovich, Silverman and Sapatinas, (1998). Contributed by Felix Abramovich and Fanis Sapatinas.
cv: See Nason, 1996.
fdr: See Abramovich and Benjamini, 1996. Contributed by Felix Abramovich.
LSuniversal: See Nason, von Sachs and Kroisandt, 1998. This is used for smoothing of a wavelet periodogram and shouldn't be used generally.
manual: specify a user supplied threshold using value to pass the value of the threshold. The value argument should be a vector. If it is of length 1 then it is replicated to be the same length as the levels vector, otherwise it is repeated as many times as is necessary to be the levels vector's length. In this way, different thresholds can be supplied for different levels. Note that the by.level option has no effect with this policy.
mannum: You decided how many of the largest (in absolute value) coefficients that you want to keep and supply this number in value.
op1: See Ogden and Parzen, 1996. Contributed by Todd Ogden.
op2: See Ogden and Parzen, 1996. Contributed by Todd Ogden.
probability: The probability policy works as follows. All coefficients that are smaller than the valueth quantile of the coefficients are set to zero. If by.level is false, then the quantile is computed for all coefficients in the levels specified by the "levels" vector; if by.level is true, then each level's quantile is estimated separately. The probability policy is pretty stupid - do not use it.
sure: See Donoho and Johnstone, 1994.
universal: See Donoho and Johnstone, 1995.

Author(s)

G P Nason

References

Various code segments detailed above were kindly donated by Felix Abramovich, Theofanis Sapatinas and Todd Ogden.

Examples

#
# Generate some test data
#
test.data <- example.1()$y
## Not run: ts.plot(test.data)
#
# Generate some noisy data
#
ynoise <- test.data + rnorm(512, sd=0.1)
#
# Plot it
#
## Not run: ts.plot(ynoise)
#
# Now take the discrete wavelet transform
# N.b. I have no idea if the default wavelets here are appropriate for
# this particular examples. 
#
ynwd <- wd(ynoise)
## Not run: plot(ynwd)
#
# Now do thresholding. We'll use a universal policy, 
# and madmad deviance estimate on the finest
# coefficients and return the threshold. We'll also get it to be verbose
# so we can watch the process.
#
ynwdT1 <- threshold(ynwd, policy="universal", dev=madmad,
		levels= nlevelsWT(ynwd)-1, return.threshold=TRUE,
	verbose=TRUE)
# threshold.wd:
# Argument checking
# Universal policy...All levels at once
# Global threshold is:  0.328410967430135 
#
# Why is this the threshold? Well in this case n=512 so sqrt(2*log(n)),
# the universal threshold,
# is equal to 3.53223. Since the noise is about 0.1 (because that's what
# we generated it to be) the threshold is about 0.353.
#
# Now let's apply this threshold to all levels in the noisy wavelet object
#
ynwdT1obj <- threshold(ynwd, policy="manual", value=ynwdT1,
	levels=0:(nlevelsWT(ynwd)-1))
#
# And let's plot it
#
## Not run: plot(ynwdT1obj)
#
# You'll see that a lot of coefficients have been set to zero, or shrunk.
#
# Let's try a Bayesian examples this time!
#
ynwdT2obj <- threshold(ynwd, policy="BayesThresh")
#
# And plot the coefficients
#
## Not run: plot(ynwdT2obj)
#
# Let us now see what the actual estimates look like
#
ywr1 <- wr(ynwdT1obj)
ywr2 <- wr(ynwdT2obj)
#
# Here's the estimate using universal thresholding
#
## Not run: ts.plot(ywr1)
#
# Here's the estimate using BayesThresh 
#
## Not run: ts.plot(ywr2)

[Package wavethresh version 4.7.2 Index]