"awmaxt"<-
function(xx, a = 0, tmax = sqrt(2 * log(length(xx))), amin = 0.01, amax = 2)
{
#  finds the marginal max likelihood estimators of a and w, given
#    data xx.   If a is set to a nonzero  value on entry then the optimisation 
#    is carried out for fixed a.  Otherwise the maximum is found over a in the range
#    (amin, amax).
#  Internally, the parametrisation is in terms of the threshold of the posterior
#    median function.  This threshold is constrained to fall in the range [0, tmax].
#
#  On exit a vector of three elements is returned, giving the optimal w, a and threshold.
#   
#
	uu <- nlminb(c(2, 0.5), marnegloglikt, gradient = marnegloglikgradt, 
		control = nlminb.control(x.tol = 0.01, rel.tol = 1e-005), lower
		 = c(0, amin), upper = c(tmax, amax), xx = xx)
	uu <- uu$parameters
	v <- c(wtfromthresh(uu[1], uu[2]), uu[2], uu[1])
	names(v) <- c("w", "a", "t")
	return(v)
}
"betafun"<-
function(x, a = 1)
{
#
#  The ratio between the standard normal and the Laplace-normal convolution
#
	x <- abs(x)
	a <- min(a, 35)
	xpa <- x + a
	xma <- x - a
	rat1 <- 1/xpa
	rat1[xpa < 35] <- pnorm( - xpa[xpa < 35])/dnorm(xpa[xpa < 35])
	rat2 <- pnorm(xma)/dnorm(xma)
	g <- (a/2) * (rat1 + rat2)
	beta1 <- 1/g
	return(beta1)
}
"betat"<-
function(w, b)
{
	sum((1 - b)/((1 - w) * b + w))
}
"cauchysmooth"<-
function(x, sig = 1, approx = T, returnthresh = F)
{
#   smooth using the Cauchy prior and the mml estimate
# of the weight
#  If approx=T, use a piecewise linear approximation to the 
#   thresholding function.
#
	x <- as.vector(x/sig)
	ww <- wchat(x)
	if(approx)
		xs <- postmedcapprox(x, w = ww)
	else xs <- sapply(x, "postmedc", w = ww)
	xsm <- xs * sig
	if(returnthresh) {
		thresh <- findcauthresh(ww)
		return(w = ww, thresh, xsm)
	}
	else return(xsm)
}
"caufun"<-
function(x, z, w)
{
# the objective function that has to be zeroed to find the 
#  posterior median.  x is the parameter, z is the data, w
# is the weight
	ez <- exp(z^2/2)
	hh <- z - x
	y <- ez * pnorm(hh) - z * ez * dnorm(hh) + (z * x - 1) * exp(z * x) * 
		pnorm( - x) - (ez - 1 + (1/w - 1) * z^2)/2
	return(y/ez)
}
"caumedfun"<-
function(x, thresh)
{
	caufun(0, thresh, x)
}
"ebaysmooth"<-
function(x, a = 0, sig = 1, returnthresh = F, posteriormean = F, hardthresh = F
	)
{
#  Find the mml estimates of the hyperparameters and then
#  find the posterior median or mean using these parameters, for 
#   a single sequence of data x
#    if a = 0, maximize over both a and w, otherwise use fixed a
#
#   if hardthresh then hard threshold with the given thresh
#
	if(a == 0) {
		pp <- awmaxt(x/sig)
		a <- pp[2]
		w <- pp[1]
		threshold <- pp[3]
	}
	else {
		w <- wthat(x/sig, a)
		if(returnthresh)
			threshold <- findexpthresh(w, a)
	}
	if(posteriormean)
		xsm <- sig * postmean(x/sig, w, a)
	else {
		if(hardthresh)
			xsm <- x * (abs(x) > sig * threshold)
		else xsm <- sig * postmed(x/sig, w, a)
	}
	if(returnthresh)
		return(w, a, threshold, xsm)
	else return(xsm)
}
"ebaywaveshrink"<-
function(x.dwt, vscale = NULL, a = 0, smooth.levels = NULL, cauchy = F, 
	postmean = F, hardthresh = F)
{
#
#  Carry out the empirical Bayes smoothing approach 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.
#
#  If cauchy=T then the cauchy prior is used
#    
#
	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)
	if(cauchy) {
		for(j in ((nlevs - slevs + 2):(nlevs + 1))) {
			x.dwt[[j]] <- cauchysmooth(as.vector(x.dwt[[j]]), sig
				 = vscale)
		}
	}
	else {
		for(j in ((nlevs - slevs + 2):(nlevs + 1))) {
			x.dwt[[j]] <- ebaysmooth(x.dwt[[j]], a, sig = vscale, 
				posteriormean = postmean, hardthresh = 
				hardthresh)
		}
	}
	x <- reconstruct(x.dwt)
	return(x)
}
"ebaywaveshrinklevdep"<-
function(x)
{
#
#  Carry out the empirical Bayes smoothing approach for the detail
#   coefficients of the wavelet transform  xdwt
#
#  use level dependent scales, chosen by mad at each scale
#
#    
#
	x.dwt <- nd.dwt(x)
	nlevs <- attributes(x.dwt)$n.levels
	for(j in (2:(nlevs + 1))) {
		x.dwt[[j]] <- ebaysmooth(x.dwt[[j]], sig = mad(x.dwt[[j]]))
	}
	x <- reconstruct(x.dwt)
	return(x)
}
"findcauthresh"<-
function(w)
{
#
#  find the threshold for weight w
#
	cauthreshfun <- function(x, ww)
	{
		caufun(0, x, ww)
	}
	uniroot(cauthreshfun, c(0.01, 6), ww = w)$root
}
"findexpthresh"<-
function(w, a = 1)
{
#
#    find the threshold that corresponds to the given weight and a.
#    
	aa <- a
	ww <- w
	if(w == 1)
		return(0)
	ft <- function(x, w, a = 1)
	{
		postmed(x, w, a) - 0.001
	}
	uniroot(ft, c(0, 10), w = ww, a = aa)$root
}
"gfun"<-
function(x, a = 1)
{
#
#  The density of the Laplace-normal convolution.
#  Take care for a greater than about 35.
#
	x <- abs(x)
	g <- (a/2) * exp(a^2/2) * (exp( - a * x) * pnorm(x - a) + exp(a * x) * 
		pnorm( - x - a))
	g[is.na(g)] <- 0
	return(g)
}
"gfungrad"<-
function(x, a = 1)
{
#
#  The gradient with respect to a of the density of the Laplace-normal convolution.
#  Take care for a greater than about 35.
#
	x <- abs(x)
	ggrad2 <- (1 + a * (a - x)) * exp(a^2/2 - a * x) * pnorm(x - a) + (1 + 
		a * (a + x)) * exp(a^2/2 + a * x) * pnorm( - x - a) - 2 * a * 
		dnorm(x)
	ggrad2[is.na(ggrad2)] <- 0
	return(ggrad2/2)
}
"marnegloglikt"<-
function(xpar, xx)
{
#  Find marginal negative log likelihood for double exponential prior
#   parametrising in terms of threshold and scale
#    xpar  is vector of parameters (thresh, a)
#    xx    is vector of data
#
	a <- xpar[2]
	w <- wtfromthresh(xpar[1], a)
	xx <- pmin(abs(xx), 200/a)
	return( - sum(log((1 - w) * dnorm(xx) + w * gfun(xx, a))))
}
"marnegloglikgradt"<-
function(xpar, xx)
{
#  Find gradient of function found by marnegloglikt
#
	a <- xpar[2]
	dd <- dwtfromthresh(xpar[1], a)
	w <- dd$w
	dw <- dd$dw
	xx <- pmin(abs(xx), 200/a)
	gf <- gfun(xx, a)
	gfg <- gfungrad(xx, a)
	dn <- dnorm(xx)
	lik <- (1 - w) * dn + w * gf
	grad1 <- dw[1] * sum((dn - gf)/lik)
	grad2 <-  - sum(((gf - dn) * dw[2] + w * gfg)/lik)
	return(c(grad1, grad2))
}
"postmean"<-
function(x, w, a = 1)
{
#
# find the posterior mean for the double exponential prior for 
#   given x, w and a, assuming the error variance
#   is 1.
#
#  only allow a < 20. 
	a <- min(a, 20)	#
#  First find the odds of zero and the shrinkage factor
#
	odzero <- ((1 - w) * betafun(x, a))/w
	wpost <- 1/(1 + odzero)	#
#  now find the posterior mean conditional on being non-zero
#
	sx <- sign(x)
	x <- abs(x)
	cp1 <- pnorm(x - a)
	dp1 <- dnorm(x - a)
	cp2 <- pnorm( - x - a)
	dp2 <- dnorm(x + a)
	ef <- exp(pmin(2 * a * x, 100))
	postmeancond <- ((x - a) * cp1 + dp1 + ef * ((x + a) * cp2 - dp2))/(cp1 +
		ef * cp2)	#
#  calculate posterior mean and return
#
	mutilde <- sx * wpost * postmeancond
	return(mutilde)
}
"postmed"<-
function(x, w, a = 1)
{
#
# find the posterior median for the double exponential prior for 
#   given x, w and a, assuming the error variance
#   is 1.
#
#  only allow a < 20. 
	a <- min(a, 20)	#
#  First find the probability of zero and of negative
#
	sx <- sign(x)
	x <- abs(x)
	odzero <- ((1 - w) * betafun(x, a))/w
	p1 <- 1/(1 + odzero)
	odpos <- (exp(-2 * a * x) * pnorm(x - a))/pnorm( - x - a)
	pneg <- p1/(1 + odpos)
	ppos <- p1 - pneg
	ppos[is.na(odpos)] <- 1	#
#  now find relevant quantile of normal, if necessary
#
	muhat <- rep(0, length(x))
	index <- (ppos > 0.5)
	if(sum(index) > 0) {
		xx <- x[index]
		qx <- pnorm(xx - a)/ppos[index]
		muhat[index] <-  - sx[index] * qnorm(qx/2, a - xx)
	}
	return(muhat)
}
"postmedc"<-
function(x, w = 0.1)
{
#
# find the posterior median function of the Cauchy prior with
#   mixing weight w
#
	ax <- abs(x)
	if(ax > 20)
		return(x - 2/x)
	flo <- caufun(0, ax, w)
	if(flo <= 0)
		return(0)
	return(sign(x) * uniroot(caufun, c(0, ax), z = ax, w = w)$root)
}
"postmedcapprox"<-
function(x, w = 0.1)
{
#
# find the approximate posterior median function of the Cauchy prior with
#   mixing weight w
#
	ax <- abs(x)
	thresh <- findcauthresh(w)
	xhat <- x - 2/x
	xhat[ax <= thresh] <- 0	#
#  carry out a linear approximation for points fairly near thresh
#
	xoff <- c(0, min(w, 0.1), 0.5, 1, 2)
	if(w > 0.4)
		xoff <- c(xoff, 3)
	if(w < 0.1)
		xoff <- sort(c(xoff, min(sqrt(w), 0.2)))
	xp <- thresh + xoff
	index <- (ax > thresh) & (ax <= max(xp))
	if(sum(index) < 5)
		xhat[index] <- sapply(ax[index], "postmedc", w = w)
	else {
		yp <- c(0, sapply(xp[-1], "postmedc", w = w))
		xhat[index] <- approx(xp, yp, ax[index])$y
	}
	xhat[index] <- sign(x[index]) * xhat[index]
	return(xhat)
}
"wchat"<-
function(xd)
{
#
#   Find the marginal maximum likelihood estimate of the mixing parameter
#    for the mixed normal prior with Cauchy tails.  It is assumed that the 
#    noise variance is equal to one.  The minimum weight is chosen to be exactly right
#    if n=1000
#
	wlo <- 0.013 * (1000/length(xd))^0.6
	whi <- 0.99
	phix <- dnorm(xd)
	gx <- xd
	index <- (xd != 0)
	gx[!index] <-  - dnorm(0)/2
	gx[index] <- (dnorm(0) - phix[index])/xd[index]^2 - phix[index]
	bcfun <- function(x, pp, gg)
	{
		sum(gg/(pp + x * gg))
	}
	flo <- bcfun(wlo, phix, gx)
	fhi <- bcfun(whi, phix, gx)
	if(flo <= 0)
		return(wlo)
	if(fhi >= 0)
		return(whi)
	uniroot(bcfun, c(wlo, whi), pp = phix, gg = gx)$root
}
"wthat"<-
function(x, a = 1)
{
#
#  Find the empirical Bayes estimate of the weight, for given data x and scale
#   parameter a, assuming that the noise standard deviation is 1
#  This is for the double exponential prior.
#
	n <- length(x)
	if(n == 1)
		return(1)
	bb <- betafun(x, a)
	wlo <- wtfromthresh(sqrt(2 * log(n)), a)
	blo <- betat(wlo, bb)
	if(blo <= 0)
		return(w = wlo)
	bhi <- betat(1, bb)
	if(bhi >= 0)
		return(w = 1)
	uu <- uniroot(betat, lower = wlo, upper = 1, f.lower = blo, f.upper = 
		bhi, b = bb)$root
	return(uu)
}
"wtfromthresh"<-
function(threshold, a = 1)
{
#
# find the weight for the mixed double exponential prior for 
#   given threshold  and scale a
#  assuming the error variance is 1.
#
#  First find the relative prob of positive and of negative
#
	x <- threshold
	ppos <- (0.5 * a * pnorm(x - a))/dnorm(x - a)
	pneg <- (0.5 * a * pnorm( - x - a))/dnorm(x + a)	#
# now solve the linear relation
#
	winv <- ppos - pneg + 1
	return(1/winv)
}
"dwtfromthresh"<-
function(threshold, a = 1)
{
#
# find the derivatives of the weight for the mixed double exponential prior for 
#   given threshold  and scale a
#  assuming the error variance is 1.
#
#  First find the relative prob of positive and of negative
#
	x <- threshold
	ppos <- pnorm(x - a)/dnorm(x - a)
	pneg <- pnorm( - x - a)/dnorm(x + a)	#
# now solve the linear relation
#
	winv <- 0.5 * a * (ppos - pneg) + 1
	w <- 1/winv	#
#  now do the same for the derivatives
#
	dwinvdt <- 0.5 * a * (2 + (x - a) * ppos - (x + a) * pneg)
	dwinvda <- 0.5 * (ppos - pneg) + 0.5 * a * ( - (x - a) * ppos - (x + a) *
		pneg)
	dwinv <- c(dwinvdt, dwinvda)
	dw <-  - w^2 * dwinv
	return(w, dw)
}