Function beta for the quasi-Cauchy prior

Description

Given a value or vector x of values, find the value(s) of the function \beta(x)=g(x)/\phi(x) - 1, where g is the convolution of the quasi-Cauchy with the normal density \phi(x).

Usage

beta.cauchy(x)

Arguments

Value

A vector the same length as x, containing the value(s) \beta(x).

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

beta.laplace

Examples

  beta.cauchy(c(-2,1,0,-4,8,50))

Function beta for the Laplace prior

Description

Given a single value or a vector of x and s, find the value(s) of the function \beta(x;s,a)=g(x;s,a)/fn(x;0,s) - 1, where fn(x;0,s) is the normal density with mean 0 and standard deviation s, and g is the convolution of the Laplace density with scale parameter a, \gamma_a(\mu), with the normal density fn(x;\mu,s) with mean mu and standard deviation s.

Usage

beta.laplace(x, s = 1, a = 0.5)

Arguments

Value

A vector of the same length as x is returned, containing the value(s) beta(x).

Note

The Laplace density is given by \gamma(u;a) = \frac{1}{2} a e^{-a|u|} and is also known as the double exponential density.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

beta.cauchy

Examples

beta.laplace(c(-2,1,0,-4,8,50), s=1)
beta.laplace(c(-2,1,0,-4,8,50), s=1:6, a=1)

Empirical Bayes thresholding on a sequence

Description

Given a sequence of data, performs Empirical Bayes thresholding, as discussed in Johnstone and Silverman (2004).

Usage

ebayesthresh(x, prior = "laplace", a = 0.5, bayesfac = FALSE, 
sdev = NA, verbose = FALSE, threshrule = "median", universalthresh = TRUE,
stabadjustment)

Arguments

Details

It is assumed that the data vector (x_1, \ldots, x_n) is such that each x_i is drawn independently from a normal distribution with mean \theta_i and variance \sigma_i^2 (\sigma_i is the same in the homogeneous case). The prior distribution of each \theta_i is a mixture with probability 1-w of zero and probability w of a given symmetric heavy-tailed distribution. The mixing weight, and possibly a scale factor in the symmetric distribution, are estimated by marginal maximum likelihood. The resulting values are used as the hyperparameters in the prior.

The true effects \theta_i can be estimated as the posterior median or the posterior mean given the data, or by hard or soft thresholding using the posterior median threshold. If hard or soft thresholding is chosen, then there is the additional choice of using the Bayes factor threshold, which is the value such thatthe posterior probability of zero is exactly half if the data value is equal to the threshold.

Value

If verbose = FALSE, a vector giving the values of the estimates of the underlying mean vector.

If verbose = TRUE, a list with the following elements:

Author(s)

Bernard Silverman

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 software for Empirical Bayes thresholding. Journal of Statistical Software, 12.

Johnstone, I. M. (2004) 'Function Estimation and Classical Normal Theory' ‘The Threshold Selection Problem’. The Wald Lectures I and II, 2004. Available from http://www-stat.stanford.edu/~imj.

Johnstone, I. M. and Silverman, B. W. (2005) Empirical Bayes selection of wavelet thresholds. Annals of Statistics, 33, 1700–1752.

The papers by Johnstone and Silverman are available from http://www.bernardsilverman.com.

See also http://www-stat.stanford.edu/~imj for further references, including the draft of a monograph by I. M. Johnstone.

See Also

tfromx, threshld

Examples

# Data with homogeneous sampling standard deviation using 
# Cauchy prior.
eb1 <- ebayesthresh(x = rnorm(100, c(rep(0,90),rep(5,10))),
                     prior = "cauchy", sdev = NA)

# Data with homogeneous sampling standard deviation using 
# Laplace prior.
eb2 <- ebayesthresh(x = rnorm(100, c(rep(0,90), rep(5,10))),
                     prior = "laplace", sdev = 1)
             
# Data with heterogeneous sampling standard deviation using 
# Laplace prior.
set.seed(123)
mu <- c(rep(0,90), rep(5,10))
sd <- c(rep(1, 40), rep(3, 60))
x  <- mu + rnorm(100, sd = sd)

# With constraints on thresholds.
eb3 <- ebayesthresh(x = x, prior = "laplace", a = NA, sdev = sd)

# Without constraints on thresholds. Observe that the estimates with
# constraints on thresholds have fewer zeroes than the estimates without
# constraints.
eb4 <- ebayesthresh(x = x, prior = "laplace", a = NA, sdev = sd,
                     universalthresh = FALSE)
print(sum(eb3 == 0))
print(sum(eb4 == 0))

# Data with heterogeneous sampling standard deviation using Laplace
# prior.
set.seed(123)
mu <- c(rep(0,90), rep(5,10))
sd <- c(rep(1, 40), rep(5,40), rep(15, 20))
x  <- mu + rnorm(100, sd = sd)

# In this example, infinity is returned as estimate when some of the
# sampling standard deviations are extremely large. However, this can
# be solved by stabilizing the data sequence before the analysis.
eb5 <- ebayesthresh(x = x, prior = "laplace", a = NA, sdev = sd)

# With stabilization.
eb6 <- ebayesthresh(x = x, prior = "laplace", a = NA, sdev = sd,
                    stabadjustment = TRUE)

Empirical Bayes thresholding on the levels of a wavelet transform.

Description

Apply an Empirical Bayes thresholding approach level by level to the detail coefficients in a wavelet transform.

Usage

ebayesthresh.wavelet(xtr, vscale = "independent", smooth.levels = Inf, 
      prior = "laplace", a = 0.5, bayesfac = FALSE, 
      threshrule = "median")

ebayesthresh.wavelet.dwt(x.dwt, vscale = "independent", smooth.levels = Inf, 
      prior = "laplace", a = 0.5, bayesfac = FALSE,
      threshrule = "median")

ebayesthresh.wavelet.wd(x.wd, vscale = "independent", smooth.levels = Inf,
      prior = "laplace", a = 0.5, bayesfac = FALSE,
      threshrule = "median")

ebayesthresh.wavelet.splus(x.dwt, vscale = "independent", smooth.levels = Inf, 
      prior = "laplace", a = 0.5, bayesfac = FALSE,
      threshrule = "median")

Arguments

Details

The routine ebayesthresh.wavelet can process a wavelet transform obtained using the routine wd in the WaveThresh R package, the routines dwt or modwt in the waveslim R package, or one of the routines (either dwt or nd.dwt) in S+WAVELETS.

Note that the wavelet transform must be calculated before the routine ebayesthresh.wavelet is called; the choice of wavelet, boundary conditions, decimated vs nondecimated wavelet, and so on, are made when the wavelet transform is calculated.

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 basic notion of processing each level of detail coefficients is easily transferred 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 = TRUE in the calls to ebayesthresh.

The main routine ebayesthresh.wavelet calls the relevant one of the routines ebayesthresh.wavelet.wd (for a transform obtained from WaveThresh), ebayesthresh.wavelet.dwt (for transforms obtained from either dwt or modwt in waveslim) or ebayesthresh.wavelet.splus (for transforms obtained from S+WAVELETS.

Value

The wavelet transform (in the same format as that supplied to the routine) of the values of the estimated regression function underlying the original data.

Author(s)

Bernard Silverman

References

Johnstone, I. M. and Silverman, B. W. (2005) Empirical Bayes selection of wavelet thresholds. Annals of Statistics, 33, 1700–1752.

See also the other references given for ebayesthresh and at http://www.bernardsilverman.com.

See Also

ebayesthresh

Weighted least squares monotone regression

Description

Given a vector of data and a vector of weights, find the monotone sequence closest to the data in the sense of weighted least squares with the given weights.

Usage

isotone(x, wt = rep(1, length(x)), increasing = FALSE)

Arguments

Details

The standard pool-adjacent-violators algorithm is used. Maximal decreasing subsequences are found within the current sequence. Each such decreasing subsequence is replaced by a constant sequence with value equal to the weighted average. Within the algorithm, the subsequence is replaced by a single point, with weight the sum of the weights within the subsequence. This process is iterated to termination. The resulting sequence is then unpacked back to the original ordering to give the weighted least squares monotone fit.

If increasing = FALSE, the original sequence is negated and the resulting estimate negated.

Value

The vector giving the best fitting monotone sequence is returned.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

wmonfromx

Posterior mean estimator

Description

Given a single value or a vector of data and sampling standard deviations (sd equals 1 for Cauchy prior), find the corresponding posterior mean estimate(s) of the underlying signal value(s).

Usage

postmean(x, s, w = 0.5, prior = "laplace", a = 0.5)
postmean.laplace(x, s = 1, w = 0.5, a = 0.5)
postmean.cauchy(x, w)

Arguments

Value

If x is a scalar, the posterior mean E(\theta|x) where \theta is the mean of the distribution from which x is drawn. If x is a vector with elements x_1, ... , x_n and s is a vector with elements s_1, ... , s_n (s_i is 1 for Cauchy prior), then the vector returned has elements E(\theta_i|x_i, s_i), where each x_i has mean \theta_i and standard deviation s_i, all with the given prior.

Note

If the quasicauchy prior is used, the argument a and s are ignored.

If prior="laplace", the routine calls postmean.laplace, which finds the posterior mean explicitly, as the product of the posterior probability that the parameter is nonzero and the posterior mean conditional on not being zero.

If prior="cauchy", the routine calls postmean.cauchy; in that case the posterior mean is found by expressing the quasi-Cauchy prior as a mixture: The mean conditional on the mixing parameter is found and is then averaged over the posterior distribution of the mixing parameter, including the atom of probability at zero variance.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

postmed

Examples

postmean(c(-2,1,0,-4,8,50), w = 0.05, prior = "cauchy")
postmean(c(-2,1,0,-4,8,50), s = 1:6, w = 0.2, prior = "laplace", a = 0.3)

Posterior median estimator

Description

Given a single value or a vector of data and sampling standard deviations (sd is 1 for Cauchy prior), find the corresponding posterior median estimate(s) of the underlying signal value(s).

Usage

postmed(x, s, w = 0.5, prior = "laplace", a = 0.5)
postmed.laplace(x, s = 1, w = 0.5, a = 0.5)
postmed.cauchy(x, w)
cauchy.medzero(x, z, w)

Arguments

Details

The routine calls the relevant one of the routines postmed.laplace or postmed.cauchy. In the Laplace case, the posterior median is found explicitly, without any need for the numerical solution of an equation. In the quasi-Cauchy case, the posterior median is found by finding the zero, component by component, of the vector function cauchy.medzero.

Value

If x is a scalar, the posterior median \mbox{med}(\theta|x) where \theta is the mean of the distribution from which x is drawn. If x is a vector with elements x_1, ... , x_n and s is a vector with elements s_1, ... , s_n (s_i is 1 for Cauchy prior), then the vector returned has elements \mbox{med}(\theta_i|x_i, s_i), where each x_i has mean \theta_i and standard deviation s_i, all with the given prior.

Note

If the quasicauchy prior is used, the argument a and s are ignored. The routine calls the approprate one of postmed.laplace or postmed.cauchy.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

postmean

Examples

postmed(c(-2,1,0,-4,8,50), w = 0.05, prior = "cauchy")
postmed(c(-2,1,0,-4,8,50), s = 1:6, w = 0.2, prior = "laplace", a = 0.3)

Find threshold from mixing weight

Description

Given a single value or a vector of weights (i.e. prior probabilities that the parameter is nonzero) and sampling standard deviations (sd equals 1 for Cauchy prior), find the corresponding threshold(s) under the specified prior.

Usage

  tfromw(w, s = 1, prior = "laplace", bayesfac = FALSE, a = 0.5)

  laplace.threshzero(x, s = 1, w = 0.5, a = 0.5)

  cauchy.threshzero(z, w)

Arguments

Details

The Bayes factor method uses a threshold such that the posterior probability of zero is exactly half if the data value is equal to the threshold. If bayesfac is set to FALSE (the default) then the threshold is that of the posterior median function given the data value.

The routine carries out a binary search over each component of an appropriate vector function, using the routine vecbinsolv.

For the posterior median threshold, the function to be zeroed is laplace.threshzero or cauchy.threshzero.

For the Bayes factor threshold, the corresponding functions are beta.laplace or beta.cauchy.

Value

The value or vector of values of the estimated threshold(s).

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

wfromx, tfromx, wandafromx

Examples

tfromw(c(0.05, 0.1), s = 1) 
tfromw(c(0.05, 0.1), prior = "cauchy", bayesfac = TRUE)

Find thresholds from data

Description

Given a vector of data and standard deviations (sd equals 1 for Cauchy prior), find the value or vector (heterogeneous sampling standard deviation with Laplace prior) of thresholds corresponding to the marginal maximum likelihood choice of weight.

Usage

  tfromx(x, s = 1, prior = "laplace", bayesfac = FALSE, a = 0.5,
         universalthresh = TRUE)

Arguments

Details

First, the routine wfromx is called to find the estimated weight. Then the routine tfromw is used to find the threshold. See the documentation for these routines for more details.

Value

The numerical value or vector of the estimated thresholds is returned.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

tfromw, wfromx

Examples

tfromx(x = rnorm(100, c(rep(0,90),rep(5,10))), prior = "cauchy")

Threshold data with hard or soft thresholding

Description

Given a data value or a vector of data, threshold the data at a specified value, using hard or soft thresholding.

Usage

threshld(x, t, hard = TRUE)

Arguments

Value

A value or vector of values the same length as x, containing the result of the relevant thresholding rule applied to x.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

ebayesthresh

Examples

threshld(-5:5, 1.4, FALSE)

Solve systems of nonlinear equations based on a monotonic function

Description

Solve a nonlinear equation or a vector of nonlinear equations based on an increasing function in a specified interval.

Usage

vecbinsolv(zf, fun, tlo, thi, nits = 30, ...)

Arguments

Details

If fun is a scalar monotone function, the routine finds a vector t the same length as zf such that, component-wise, fun(t) = zf, where this is possible within the interval \code{(tlo,thi)}. The relevant value returned is the nearer extreme to the solution if there is no solution in the specified range for any particular component of zf. The routine will also work if fun is a vector of monotone functions, allowing different functions to be considered for different components. The interval over which the search is conducted has to be the same for each component.

The accuracy of the solution is determined by the number of binary subdivisions; if nits=30 then the solution(s) will be accurate to about 9 orders of magnitude less than the length of the original interval (tlo, thi).

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

Find weight and scale factor from data if Laplace prior is used.

Description

Given a vector of data and a single value or vector of sampling standard deviations, find the marginal maximum likelihood choice of both weight and scale factor under the Laplace prior.

Usage

wandafromx(x, s = 1, universalthresh = TRUE)
negloglik.laplace(xpar, xx, ss, tlo, thi)

Arguments

Details

The parameters are found by marginal maximum likelihood.

The search is over weights corresponding to threshold t_i in the range [0, s_i \sqrt{2 \log n}] if universalthresh=TRUE, where n is the length of the data vector and (s_1, ... , s_n) is the vector of sampling standard deviation of data (x_1, ... , x_n); otherwise, the search is over [0,1].

The search uses a nonlinear optimization routine (optim in R) to minimize the negative log likelihood function negloglik.laplace. The range over which the scale factor is searched is (0.04, 3). For reasons of numerical stability within the optimization, the prior is parametrized internally by the threshold and the scale parameter.

Value

A list with values:

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

wfromx, tfromw

Examples

wandafromx(rnorm(100, c(rep(0,90),rep(5,10))), s = 1)

Mixing weight from posterior median threshold

Description

Given a value or vector of thresholds and sampling standard deviations (sd equals 1 for Cauchy prior), find the mixing weight for which this is(these are) the threshold(s) of the posterior median estimator. If a vector of threshold values is provided, the vector of corresponding weights is returned.

Usage

wfromt(tt, s = 1, prior = "laplace", a = 0.5)

Arguments

Value

The numerical value or vector of values of the corresponding weight is returned.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

tfromw

Examples

wfromt(c(2,3,5), prior = "cauchy" )

Find Empirical Bayes weight from data

Description

Suppose the vector (x_1, \ldots, x_n) is such that x_i is drawn independently from a normal distribution with mean \theta_i and standard deviation s_i (s_i equals 1 for Cauchy prior). The prior distribution of the \theta_i is a mixture with probability 1-w of zero and probability w of a given symmetric heavy-tailed distribution. This routine finds the marginal maximum likelihood estimate of the parameter w.

Usage

wfromx(x, s = 1, prior = "laplace", a = 0.5, universalthresh = TRUE)

Arguments

Details

The weight is found by marginal maximum likelihood.

The search is over weights corresponding to threshold t_i in the range [0, s_i \sqrt{2 \log n}] if universalthresh=TRUE, where n is the length of the data vector and (s_1, ... , s_n) (s_i is 1 for Cauchy prior) is the vector of sampling standard deviation of data (x_1, ... , x_n); otherwise, the search is over [0, 1].

The search is by binary search for a solution to the equation S(w)=0, where S is the derivative of the log likelihood. The binary search is on a logarithmic scale in w.

If the Laplace prior is used, the scale parameter is fixed at the value given for a, and defaults to 0.5 if no value is provided. To estimate a as well as w by marginal maximum likelihood, use the routine wandafromx.

Value

The numerical value of the estimated weight.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

wandafromx, tfromx, tfromw, wfromt

Examples

wfromx(x = rnorm(100, s = c(rep(0,90),rep(5,10))), prior = "cauchy")

Find monotone Empirical Bayes weights from data.

Description

Given a vector of data, find the marginal maximum likelihood choice of weight sequence subject to the constraints that the weights are monotone decreasing.

Usage

wmonfromx(xd, prior = "laplace", a = 0.5, tol = 1e-08, maxits = 20)

Arguments

Details

The weights is found by marginal maximum likelihood. The search is over weights corresponding to thresholds in the range [0, \sqrt{2 \log n}], where n is the length of the data vector.

An iterated least squares monotone regression algorithm is used to maximize the log likelihood. The weighted least squares monotone regression routine isotone is used.

To turn the weights into thresholds, use the routine tfromw; to process the data with these thresholds, use the routine threshld.

Value

The vector of estimated weights is returned.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

wfromx, isotone

Estimation of a parameter in the prior weight sequence in the EbayesThresh paradigm.

Description

Suppose a sequence of data has underlying mean vector with elements \theta_i. Given the sequence of data, and a vector of scale factors cs and a lower limit pilo, this routine finds the marginal maximum likelihood estimate of the parameter zeta such that the prior probability of \theta_i being nonzero is of the form median(pilo, zeta*cs, 1).

Usage

zetafromx(xd, cs, pilo = NA, prior = "laplace", a = 0.5)

Arguments

Details

An exact algorithm is used, based on splitting the range up for zeta into subintervals over which no element of zeta*cs crosses either pilo or 1.

Within each of these subintervals, the log likelihood is concave and its maximum can be found to arbitrary accuracy; first the derivatives at each end of the interval are checked to see if there is an internal maximum at all, and if there is this can be found by a binary search for a zero of the derivative.

Finally, the maximum of all the local maxima over these subintervals is found.

Value

A list with the following elements:

Note

Once the maximizing zeta and corresponding weights have been found, the thresholds can be found using the program tfromw, and these can be used to process the original data using the routine threshld.

Author(s)

Bernard Silverman

References

See ebayesthresh and http://www.bernardsilverman.com

See Also

tfromw, threshld, wmonfromx, wfromx