treethresh (version 0.1-11)

treethresh: Compute optimal thresholding partition for a sequence, matrix, or array of data.

Description

This function carries out the tree-based thresholding algorithm described in section 3 of Evers and Heaton (2009).

Usage

treethresh(data, beta, criterion="score", control=list(),
           rho=sys.frame(sys.parent()))

Arguments

data

An array (or an object coercible to an array, i.e. a vector or matrix) containing the data. The data is assumed to have noise of unit variance, thus the data needs to be rescaled a priori (e.g. in the case of wavelet coefficients using function estimate.sdev).

beta

Instead of using the original data, one can call wtthresh using the \(\beta_i\) instead of the observed data. These can be computed using beta.laplace.

criterion

The criterion to be used. Can be "score" (default) for using the score test, "likelihood" for using the likelihood ratio test (slower), "heuristic" for using a heuristic criterion based on the original data, or a user-specified function that computes the goodness of a split. This function should take four arguments (which should be self-explanatory), left_data, left_betas, right_data, and right_betas.

control

A list that allows the user to tweak the behaviour of treethresh. It can contain the following elements:

max.depth

The maximum depth of the tree. Defaults to 10.

minimum.width

The minimum width of a region of the partitions. This setting avoids creating too small regions. Defaults to 3.

minimum.size

The minimum size of a region of the partitions. This setting avoids creating too small regions. Defaults to \(5^d\), where \(d\) is the dimension of the arras.

lr.signif

If the p-value of the corresponding likelihood ratio test is larger than 1-lr.signif a split will be discarded. Defaults to 0.5.

absolute.improvement

The minimum absolute improvement of the above criterion necessary such that a split is retained. Defaults to -Inf, i.e. deactivated.

relative.improvement

The minimum relative improvement of the above criterion necessary such that a split is retained. Defaults to -Inf, i.e. deactivated.

absolute.criterion

The minimum value of the above criterion necessary such that a split is retained. Defaults to 0, i.e. deactivated.

a

The parameter \(a\) of the Laplace distribution \(\gamma(\mu)\propto \exp(-a|\mu|) \) corresponding to the signal. Defaults to 0.5.

beta.max

The maximum value of \(\beta\). Defaults to 1e5.

max.iter

The maximum number of iterations when computing the estimate of the weight \(w\) in a certain region. Defaults to 30.

tolerance.grad

The estimate of the weight \(w\) in a certain region is considered having converged, if the gradient of the likelihood is less than tolerance.grad. Defaults to 1e-8.

tolerance

The estimate of the weight \(w\) in a certain region is considered having converged, if the estimates of the weight \(w\) change less than tolerance. Defaults to 1e-6.

rho

The environment used to evaluate the user-speficied criterion function if one is supplied). (You want to change this argument only in very rare circumstances).

Value

treethresh returns an object of the class c("treethresh"), which is a list containing the following elements:

splits

A table describing the detailed structure of the fitted tree together with the local loglikelihoods required for the pruning.

membership

An array of the same dimension as data or beta indicating to which region each entry of the array of data belongs.

beta

The values of \(\beta\) for each observation / coefficient.

data

The data used.

criterion

The criterion used to decide on splits (see argument criterion).

control

The control list of tuning options used (see argument control).

References

Evers, L. and Heaton T. (2009) Locally Adaptive Tree-Based Thresholding, Journal of Computational and Graphical Statistics 18 (4), 961-977. Evers, L. and Heaton T. (2017) Locally Adaptive Tree-Based Thresholding, Journal of Statistical Software, Code Snippets, 78(2), 1-22.

Examples

# (1) Create a vector with the probabilities of a signal being present
w.true <- c(rep(0.1,400),rep(0.7,300),rep(0.1,300))

# (2) Generate the signal
mu <- numeric(length(w.true))
non.zero.entry <- runif(length(mu))<w.true
num.non.zero.entries <- sum(non.zero.entry)
mu[non.zero.entry] <- rexp(num.non.zero.entries,rate=0.5)*
                         sample(c(-1,1),num.non.zero.entries,replace=TRUE)

# (3) Split graphics device
par(mfrow=c(2,2))

# (3) Draw the true signal (signal present in red)
plot(mu,col=non.zero.entry+1)
title("True signal")

# (4) Add noise to the signal
x <- mu + rnorm(length(mu))

# (5) Plot the noisy signal (signal present in red)
plot(x,col=non.zero.entry+1)
title("Noisy signal")

# (6) Carry out the tree-based thresholding
tt <- treethresh(x)

# (7) Prune the tree
tt.pruned <- prune(tt)

# (8) Threshold the signal according to the pruned tree
mu.hat <- thresh(tt.pruned)

# (9) Print the denoised signal
plot(mu.hat,col=non.zero.entry+1)
title("Denoised signal")

# (10) Add solid lines for splits (lines removed by the pruing are dashed)
abline(v=tt$split[,"pos"],lty=2)
abline(v=tt.pruned$split[,"pos"],lty=1)