"simulation1"<-
function(nreps = 100)
{
	kvals <- c(5, 50, 500)
	muvals <- c(3, 4, 5, 7)
	results <- array(NA, dim = c(11, 4, 3), dimnames = list(c("exponential",
		"cauchy", "postmean", "exphard", "SURE", "adapt", "FDR q=0.01", 
		"FDR q=0.1", "FDR q=0.4", "universal soft", "universal hard"), 
		as.character(muvals), as.character(kvals)))	
	#  simulate matrix of errors
#
	set.seed(153)
	errmat <- matrix(rnorm(nreps * 1000), nrow = 1000)	#
#  now conduct individual simulations
#
	for(jk in (1:3)) {
		for(jmu in (1:4)) {
			zmean <- rep(0, 1000)
			zmean[1:kvals[jk]] <- muvals[jmu]
			z <- matrix(rep(zmean, nreps), nrow = 1000) + errmat
			results[1, jmu, jk] <- simcomp1("ebaysmooth", zmean, z)
			results[2, jmu, jk] <- simcomp1("cauchysmooth", zmean, 
				z)
			results[3, jmu, jk] <- simcomp1("ebaysmooth", zmean, z, 
				posteriormean = T)
			results[4, jmu, jk] <- simcomp1("ebaysmooth", zmean, z, 
				hardthresh = T)
			results[5, jmu, jk] <- simcomp1("suresmooth", zmean, z, 
				adapt = F)
			results[6, jmu, jk] <- simcomp1("suresmooth", zmean, z, 
				adapt = T)
			results[7, jmu, jk] <- simcomp1("fdrsmooth", zmean, z, 
				q = 0.01)
			results[8, jmu, jk] <- simcomp1("fdrsmooth", zmean, z, 
				q = 0.1)
			results[9, jmu, jk] <- simcomp1("fdrsmooth", zmean, z, 
				q = 0.4)
			results[10, jmu, jk] <- simcomp1("universalsmooth", 
				zmean, z, soft = T)
			results[11, jmu, jk] <- simcomp1("universalsmooth", 
				zmean, z, soft = F)
		}
	}
	return(results)
}
"simulation1a"<-
function(nreps = 100, avals = c(1, 0.5, 0.2, 0.1, 0))
{
#  Compare results for optimisation with fixed a to more general optimisation
#
	navals <- length(avals)
	kvals <- c(5, 50, 500)
	muvals <- c(3, 4, 5, 7)
	results <- array(NA, dim = c(navals, 4, 3), dimnames = list(
		as.character(avals), as.character(muvals), as.character(kvals))
		)	#  simulate matrix of errors
#
	set.seed(153)
	errmat <- matrix(rnorm(nreps * 1000), nrow = 1000)	#
#  now conduct individual simulations
#
	for(jk in (1:3)) {
		for(jmu in (1:4)) {
			zmean <- rep(0, 1000)
			zmean[1:kvals[jk]] <- muvals[jmu]
			z <- matrix(rep(zmean, nreps), nrow = 1000) + errmat
			for(ja in (1:navals)) {
				results[ja, jmu, jk] <- simcomp1("ebaysmooth", 
				  zmean, z, a = avals[ja])
			}
		}
	}
	return(results)
}
"simulation2"<-
function(nreps = 100, nondecimated = T)
{
	results <- array(NA, dim = c(9, 4, 2, nreps))
	signals <- c("bumps", "blocks", "doppler", "heavisine")
	methods <- c("exp", "cau", "postmean", "exphard", "sure", "sure6", 
		"univ", "spline", "tukey")
	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,  ] <- simstudy(signals[i], snrs[j], 
				nreps, nondecimated)
			cat("+")
		}
	}
	return(results)
}
"simulation3"<-
function(nreps = 100)
{
#
#  Carries out simulation on the same data as simulation2, using the Cai/Silverman
#  NeighBlock and NeighCoeff approaches
#
	results <- array(NA, dim = c(2, 4, 2, nreps))
	signals <- c("bumps", "blocks", "doppler", "heavisine")
	methods <- c("NeighBlock", "NeighCoeff")
	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,  ] <- blocksimstudy(signals[i], snrs[j],
				nreps, nondecimated)
			cat("+")
		}
	}
	return(results)
}
"table2"<-
function(zs = sim2.res, tex = F)
{
#   takes results of simulation2() or simulation3() and
#   produces Table 2.   Can either produce
#   an array or else latex output
	zr <- apply(zs, c(1, 2, 3), mean)
	if(!tex) {
		return(zr)
	}
	zt <- round(cbind(zr[,  , 1], zr[,  , 2]))
	zt <- apply(zt, 1, paste, collapse = " & ")
	zt <- paste(names(zt), zt, sep = " & ", collapse = " \\  ")
	return(zt)
}
"table3"<-
function(zs = sim2.res, tex = F)
{
#   takes results of simulation2() and
#   produces Table 3.   Can either produce
#   an array or else latex output
	zr <- zs[2:7,  ,  ,  ]
	for(j in 1:6)
		zr[j,  ,  ,  ] <- zr[j,  ,  ,  ] - zs[1,  ,  ,  ]
	zrdif <- apply(zr, c(1, 2, 3), mean)
	zrse <- sqrt(apply(zr, c(1, 2, 3), var)/(dim(zr)[4]))
	zr <- zrdif/zrse
	if(!tex) {
		return(zr)
	}
	zr <- cbind(zr[,  , 1], zr[,  , 2])
	index <- (abs(zr) <= 4)
	zr[index] <- 10 * zr[index]
	zt <- round(zr)
	zt[index] <- zt[index]/10
	zt <- apply(zt, 1, paste, collapse = " & ")
	zt <- paste(names(zt), zt, sep = " & ", collapse = " \\  ")
	return(zt)
}
"simcomp"<-
function(zdata, ztrue, nondecimated = T, ...)
{
#  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, wavelet shrinkage using ebayes and the hard threshold chosen by
#   that method, sure using two primary resolution levels,  
#   universal thresholding, spline smoothing and the default Splus smooth.  
#
	nn <- length(zdata)
	error <- rep(NA, 9)
	if(nondecimated)
		z.dwt <- nd.dwt(zdata, ...)
	else z.dwt <- dwt(zdata, ...)
	nlevs <- attributes(z.dwt)$n.levels
	sigest <- mad(z.dwt[[nlevs + 1]])
	error[1] <- vecnorm(ebaywaveshrink(z.dwt, vscale = sigest) - ztrue)
	error[2] <- vecnorm(ebaywaveshrink(z.dwt, vscale = sigest, cauchy = T) - 
		ztrue)
	error[3] <- vecnorm(ebaywaveshrink(z.dwt, vscale = sigest, postmean = T
		) - ztrue)
	error[4] <- vecnorm(ebaywaveshrink(z.dwt, vscale = sigest, hardthresh
		 = T) - ztrue)
	error[5] <- vecnorm(waveshrink(z.dwt, vscale = sigest, shrink.rule = 
		"sure") - ztrue)
	error[6] <- vecnorm(waveshrink(z.dwt, vscale = sigest, shrink.rule = 
		"sure", smooth.levels = 6) - ztrue)
	error[7] <- vecnorm(waveshrink(z.dwt, vscale = sigest, shrink.rule = 
		"universal") - ztrue)
	error[8] <- vecnorm(smooth.spline(zdata)$y - ztrue)
	error[9] <- vecnorm(as.vector(smooth(zdata)) - ztrue)
	return(error^2)
}
"simcomp1"<-
function(method, zmean, z, ...)
{
	nreps <- dim(z)[2]
	zs <- apply(z, 2, method, ...)
	mse <- vecnorm(zs - zmean)^2/nreps
	cat(method, "\n")
	return(mse)
}
"simstudy"<-
function(signal = "bumps", snr = 3, nreps = 100, nondecimated = T)
{
	set.seed(153)
	ztrue <- make.signal(signal)
	ztrue <- (snr * ztrue)/sqrt(var(ztrue))
	errs <- matrix(NA, nrow = 9, ncol = nreps)
	for(j in (1:nreps)) {
		zdata <- rnorm(1024, ztrue)
		errs[, j] <- simcomp(zdata, ztrue, nondecimated)
		cat(".")
	}
	return(errs)
}
"blocksimstudy"<-
function(signal = "bumps", snr = 3, nreps = 100, nondecimated = T)
{
	set.seed(153)
	ztrue <- make.signal(signal)
	ztrue <- (snr * ztrue)/sqrt(var(ztrue))
	errs <- matrix(NA, nrow = 2, ncol = nreps)
	for(j in (1:nreps)) {
		zdata <- rnorm(1024, ztrue)
		errs[1, j] <- (vecnorm(NeighBlock(zdata) - ztrue))^2
		errs[2, j] <- (vecnorm(NeighCoeff(zdata) - ztrue))^2
		cat(".")
	}
	return(errs)
}
"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)
}
"suresmooth"<-
function(x, adapt = T, returnthresh = F)
{
#
#  find the SURE threshold and apply it to the sequence x
#  also find carry out the ADAPT smoothing if adapt=T
#
	n <- length(x)	#
#  first do the adapt test 
#
	if(adapt & (sum(x^2)/n <= 1 + (log(n)/log(2))^1.5/sqrt(n)))
		thresh <- sqrt(2 * log(n))
	else {
#
#  the sure calculation
#
		sx <- sort(x^2)
		sureseq <- n - 2 * (1:n) + seq(from = n - 1, to = 0, by = -1) * 
			sx + cumsum(sx)
		thresh <- sqrt(min(sx[sureseq == min(sureseq)], 2 * log(n)))
	}
	xhat <- sign(x) * pmax(0, abs(x) - thresh)
	if(returnthresh)
		return(thresh)
	else return(xhat)
}
"universalsmooth"<-
function(x, soft = T)
{
	n <- length(x)
	tt <- sqrt(2 * log(n))
	sx <- sign(x)
	ax <- abs(x)
	if(soft)
		xt <- sx * pmax(0, ax - tt)
	else {
		xt <- x
		xt[ax <= tt] <- 0
	}
	return(xt)
}
"NeighBlock"<-
function(y, wavelet = "s8", j0 = NULL, extension = NULL)
{
	n <- length(y)
	J <- ceiling(log(n)/log(2))
	L0 <- floor(log(2^J)/2)
	L1 <- max(1, floor(L0/2))
	L <- L0 + 2 * L1
	if(is.null(j0))
		j0 <- ceiling(log(log(2^J))/log(2)) + 1
	j0 <- max(j0, floor(log(log(2^J))/log(2)))	
	#Use MAD to estimate the noise level of the signal. 
	noise.level <- median(abs(diff(y) - median(diff(y))))/(sqrt(2) * 0.6745
		)
	thresh <- 4.505241 * L * noise.level^2	#check the boundary.
	nlevels <- J - j0
	if(2^J > n) {
		add.left <- floor((2^J - n)/2)
		add.right <- ceiling((2^J - n)/2)
		if(is.null(extension))
			extension <- "reflection"
		if(extension == "reflection") {
			if(add.left == 0)
				y <- c(y, y[n:(n - add.right + 1)])
			else y <- c(y[add.left:1], y, y[n:(n - add.right + 1)])
		}
		else if(extension == "periodic") {
			if(add.left == 0)
				y <- c(y, y[1:add.right])
			else y <- c(y[(n - add.left + 1):n], y, y[1:add.right])
		}
		else if(extension == "zero") {
			add.left <- 0
			add.right <- 2^J - n
			y <- c(y, rep(0, add.right))
		}
		else stop("Available extension rule: periodic, reflection, zero"
				)
	}
	ydwt <- dwt(y, wavelet, n.levels = nlevels)
	for(level in 1:nlevels) {
		y <- as.vector(ydwt[[paste("d", level, sep = "")]])
		m <- length(y)
		nb <- floor(m/L0)
		yy <- c(y[(m - L1 + 1):m], y, y[1:L1])
		SS <- rep(0, nb)
		for(b in 1:nb) {
			SS[b] <- sum(yy[((b - 1) * L0 + 1):((b - 1) * L0 + L)]^
				2)
		}
		factor <- pmax(0, 1 - thresh/SS)
		factor <- rep(factor, rep(L0, nb))
		y[1:(nb * L0)] <- y[1:(nb * L0)] * factor
		if(m > nb * L0) {
			Lb <- floor((L - (m - nb * L0))/2)
			SSL <- sum(c(y[(m - L + Lb + 1):m], y[1:Lb])^2)
			y[(nb * L0 + 1):m] <- y[(nb * L0 + 1):m] * max(0, 1 - 
				thresh/SSL)
		}
		ydwt[[paste("d", level, sep = "")]] <- y
	}
	if(2^J > n)
		output <- reconstruct(ydwt)[(add.left + 1):(add.left + n)]
	else output <- reconstruct(ydwt)
	output
}
"NeighCoeff"<-
function(y, wavelet = "s8", shrink.rule = "soft", j0 = NULL, extension = NULL)
{
	n <- length(y)
	J <- ceiling(log(n)/log(2))
	if(is.null(j0))
		j0 <- ceiling(log(log(2^J))/log(2)) + 1
	j0 <- max(j0, 2)	#Use MAD to estimate the noise level of the signal. 
	noise.level <- median(abs(diff(y) - median(diff(y))))/(sqrt(2) * 0.6745
		)
	thresh <- 2 * log(n) * noise.level^2
	nlevels <- J - j0
	add <- 2^J - n
	if(add > 0) {
		add.left <- floor(add/2)
		add.right <- add - add.left
		if(is.null(extension))
			extension <- "reflection"
		if(extension == "reflection") {
			if(add.left == 0)
				y <- c(y, y[n:(n - add.right + 1)])
			else y <- c(y[add.left:1], y, y[n:(n - add.right + 1)])
		}
		else if(extension == "periodic") {
			if(add.left == 0)
				y <- c(y, y[1:add.right])
			else y <- c(y[(n - add.right + 1):n], y, y[1:add.left])
		}
		else if(extension == "zero") {
			add.left <- 0
			add.right <- 2^J - n
			y <- c(y, rep(0, add.right))
		}
		else stop("Available extension rule: periodic, reflection, zero"
				)
	}
	ydwt <- dwt(y, wavelet, n.levels = nlevels)
	for(level in 1:nlevels) {
		y <- as.vector(ydwt[[paste("d", level, sep = "")]])
		m <- length(y)
		yy <- c(y[m], y, y[1])
		SS <- rep(0, m)
		for(b in 1:m) {
			SS[b] <- sum((yy[b:(b + 2)])^2)
		}
		if(shrink.rule == "soft") {
			ydwt[[paste("d", level, sep = "")]] <- y * pmax(0, 1 - 
				thresh/SS)
		}
		else if(shrink.rule == "hard") {
			y[SS <= thresh] <- 0
			ydwt[[paste("d", level, sep = "")]] <- y
		}
		else stop("Available shrink.rule: soft, hard")
	}
	if(2^J > n)
		output <- reconstruct(ydwt)[(add.left + 1):(add.left + n)]
	else output <- reconstruct(ydwt)
	output
}