"fdrsimcomp"<-
function(zdata, ztrue, nondecimated = T, qvec = c(0.01, 0.1, 0.4), ...)
{
#  perform a comparison between various methods of estimating ztrue
#   methods are :  wavelet shrinkage using ebayes, wavelet shrinkage
#   using the cauchy prior, wavelet shrinkage using ebayes and the posterior
#   mean, sure using two primary resolution levels,  
#   universal thresholding, spline smoothing and the default Splus smooth.  
#
	nq <- length(qvec)
	error <- rep(NA, nq)
	if(nondecimated)
		z.dwt <- nd.dwt(zdata, ...)
	else z.dwt <- dwt(zdata, ...)
	for(jq in (1:nq)) {
		error[jq] <- vecnorm(fdrwaveshrink(z.dwt, q = qvec[jq]) - ztrue
			)
	}
	return(error^2)
}
"fdrsimstudy"<-
function(signal = "bumps", snr = 3, nreps = 100, qvec = c(0.01, 0.1, 0.4), 
	nondecimated = T)
{
	set.seed(153)
	ztrue <- make.signal(signal)
	ztrue <- (snr * ztrue)/sqrt(var(ztrue))
	errs <- matrix(NA, nrow = length(qvec), ncol = nreps)
	dimnames(errs) <- list(as.character(qvec), NULL)
	for(j in (1:nreps)) {
		zdata <- rnorm(1024, ztrue)
		errs[, j] <- fdrsimcomp(zdata, ztrue, nondecimated, qvec)
		cat(".")
	}
	return(errs)
}
"fdrsimulation"<-
function(nreps = 100, nondecimated = T)
{
	results <- array(NA, dim = c(4, 4, 2, nreps))
	signals <- c("bumps", "blocks", "doppler", "heavisine")
	methods <- c("FDR0.01", "FDR0.05", "FDR0.1", "FDR0.4")
	snrnames <- c("High", "Low")
	snrs <- c(3, 7)
	dimnames(results) <- list(methods, signals, snrnames, NULL)
	for(i in (1:4)) {
		for(j in (1:2)) {
			results[, i, j,  ] <- fdrsimstudy(signals[i], snrs[j], 
				nreps, qvec = c(0.01, 0.05, 0.1, 0.4), 
				nondecimated)
			cat("+")
		}
	}
	return(results)
}
"fdrsmooth"<-
function(x, q = 0.05)
{
#  implement the FDR procedure
#
	n <- length(x)
	y <- rev(sort(abs(x)))
	t <-  - qnorm(((1:n) * q)/(2 * n))
	tfdr <- min(t[y >= t], t[1])
	xfdr <- x
	xfdr[abs(x) <= tfdr] <- 0
	return(xfdr)
}
"fdrwaveshrink"<-
function(x.dwt, q = 0.05, vscale = NULL, smooth.levels = NULL)
{
#
#  Carry out the FDR thresholding method for the detail
#   coefficients of the wavelet transform  xdwt
#
#  If smooth.levels is set to a smaller number than n.levels then only smooth.levels
#  levels are processed.
#
#
	nlevs <- attributes(x.dwt)$n.levels
	if(is.null(smooth.levels))
		smooth.levels <- nlevs
	if(is.null(vscale))
		vscale <- mad(x.dwt[[nlevs + 1]])
	slevs <- min(nlevs, smooth.levels)
	for(j in ((nlevs - slevs + 2):(nlevs + 1))) {
		x.dwt[[j]] <- fdrsmooth(x.dwt[[j]], q)
	}
	x <- reconstruct(x.dwt)
	return(x)
}