ebayesthresh_wavelet

Empirical Bayes Wavelet Thresholding

DESCRIPTION

Performs Empirical Bayes Thresholding on the levels of a wavelet transform

USAGE

x = ebayesthresh_wavelet(xdwt,L,qmf,vscale,smooth_levels,prior,a,bayesfac,isotone,threshrule)

REQUIRED ARGUMENTS

xdwt	wavelet coefficients

OPTIONAL ARGUMENTS

L	coarse approximation level for dwt (see Wavelab Toolbox)
qmf	quadrature mirror filter (see Wavelab Toolbox)
vscale	Controls the scale used at different levels of the transform. If vscale is a scalar quantity, then it will be assumed that the wavelet coefficients at every level have this standard deviation. If vscale = 'independent', the standard deviation will be estimated from the highest level of the wavelet transform and will then be used for all levels processed. If vscale='level', then the standard deviation will be estimated separately for each level processed, allowing standard deviation that is level-dependent.
smooth_levels	the number of levels to be processed, if less than the number of levels of detail calculated by the wavelet transform
prior	the prior to be used conditional on the mean being nonzero. Possible values are 'laplace' for the Laplace prior and 'cauchy' for the quasi-Cauchy prior.
a	the scale factor if the Laplace prior is used. If, on entry, a=[] and prior='laplace', then the scale parameter a will be estimated by marginal maximum likelihood, separately for each level of the transform. If the quasi-Cauchy prior is used then this argument is ignored.
bayesfac	if bayesfac=1 then whenever a threshold is explicitly calculated, the Bayes factor threshold will be used, viz the value of the data such that the posterior probability that the underlying mean is zero is exactly 0.5. Otherwise the posterior median threshold will be used.
isotone	0 or 1 (switch for monotone marginal maximum likelihood estimation). Default = 0.
threshrule	specifies the thresholding rule to be applied to the data. Possible values are 'median' (use the posterior median); 'mean' (use the posterior mean); 'hard' (carry out hard thresholding); 'soft' (carry out soft thresholding).

VALUE

x A vector giving the values of the estimated regression function underlying the original data

BACKGROUND

Given a wavelet transform obtained using the routine (FWT_PO) in WAVELAB, apply an Empirical Bayes thresholding approach to the detail coefficients in the transform, as discussed in Johnstone and Silverman (2002a,b). Only the standard wavelet transform can be processed. After the thresholding, the wavelet transform is inverted to obtain the reconstruction of the estimated regression function underlying the original data. Apart from some housekeeping to estimate the standard deviation if necessary, and to determine the number of levels to be processed, the main part of the routine is a call, for each level, to the smoothing routine ebayesthresh. The final step is a call to the WAVELAB routine IWT_PO to invert the transform. The basic notion of processing each level of detail coefficients is easily transfered to transforms constructed using other wavelet software. Similarly, it is straightforward to modify the routine to give other details of the wavelet transform, if necessary using the option verbose = 1 in the calls to ebayesthresh.

REFERENCES

Johnstone, I. M. and Silverman, B. W. (2004) Needles and straw in haystacks: Empirical Bayes estimates of possibly sparse sequences. Annals of Statistics, 32, 1594–1649.
Johnstone, I. M. and Silverman, B. W. (2004) EbayesThresh: R and S-PLUS software for Empirical Bayes thresholding. Submitted for publication.
Johnstone, I. M. and Silverman, B. W. (2005) Empirical Bayes selection of wavelet thresholds. Annals of Statistics, 33?, to appear.
The papers by Johnstone and Silverman are available from http://www.bernardsilverman.com.

SEE ALSO

ebayesthresh

EXAMPLES

n=256;
z=MakeSignal('HeaviSine',n);
y=z + 0.8*randn(size(z));
qmf = MakeONFilter('Symmlet',8);
L=3;
ydwt=FWT_PO(y,L,qmf);
yr=ebayesthresh_wavelet(ydwt,3,qmf,'independent',3,'laplace',[],0,0,'median');
plot((1:n)/n,yr,'r',(1:n)/n,z,'black',(1:n)/n,y,'o')