The Γ-minimax paradigm, originally proposed by Robbins [29], deals with the problem of selecting decision rules in tasks of statistical inference. The Γ-minimax approach falls between the Bayes paradigm, which selects procedures that work well “on average a posteriori” for a specific prior, and the minimax paradigm, which guards against least favorable outcomes, however unlikely. This approach has evolved from seminal papers in the fifties [29, 21] and early sixties, through an extensive research on foundations and parametric families in the seventies, to a branch of Bayesian robustness theory, in the eighties and nineties. A comprehensive discussion of the Γ-minimax approach can be found in [5, 6]. Also see [8], a fundamental contribution.

The Γ-minimax paradigm incorporates the prior information about the statistical model by a family of plausible priors, denoted by Γ, rather than by a single prior. Elicitation of “prior families” is often encountered in practice. Given the family of priors, the decision maker selects a rule that is optimal with respect to the least favorable prior in the family.

We note that least favorable priors typically do not exist in unbounded parameter spaces. For example, there is no least favorable prior for estimation of unconstrained normal or Poisson means. There is a least favorable prior for estimating a binomial parameter, and this is a compact parameter space. But no amount of intuition can predict the least favorable prior for the binomial. And, as soon as the loss is changed to a normalized squared error, the least favorable prior changes. There are no links or intuition regarding least favorable priors that are portable. It is too delicate and problem specific. For example, in a quite general setup, Clarke and Barron in [13] show that Jeffreys’ prior is asymptotically least favorable when the loss is Kullback-Leibler distance, thus providing theoretical justification for reference priors.

Inference of this kind is often interpreted in terms of game theory. Formally, if D is a set of all decision rules and Γ is a family of prior distributions over the parameter space Θ, then a rule δ ∗ ∈ D is Γ-minimax if (1.1) inf δ ∈ D sup π ∈ Γ r ( π , δ ) = sup π ∈ Γ r ( π , δ ∗ ) = r ( π ∗ , δ ∗ ) , where r ( π , δ ) = E θ E X | θ L ( θ , δ ) = E θ R ( θ , δ ) is the Bayes risk under the loss L ( θ , δ ). Here R ( θ , δ ) = E θ X | θ L ( θ , δ ) denotes the frequentist risk of δ, π ∗ is the least favorable prior, and L ( θ , δ ) is the loss function, often the squared error loss, ( θ − δ ) 2 . Note that when Γ is the set of all priors, the Γ-minimax rule coincides with the minimax rule, if a minimax rule exists; when Γ contains a single prior, then the Γ-minimax rule coincides with Bayes’ rule with respect to that prior. When the decision problem, viewed as a statistical game, has a value, that is, when inf δ ∈ D sup π ∈ Γ ≡ sup π ∈ Γ inf δ ∈ D , then the Γ-minimax solution coincides with the Bayes rule with respect to the least favorable prior. For the interplay between the Γ-minimax and Bayesian paradigms, see [6]. A review on Γ-minimax estimation can be found in [32].

In a nutshell, the Γ-minimax is a philosophical compromise between subjective Bayes with a single prior and the conservative approach of not having a prior at all.

We apply the Γ-minimax approach to the classical nonparametric regression problem (1.2) y i = f ( t i ) + σ ε i , i = 1 , … , n , where t i , i = 1 , … , n, are deterministic equispaced design points on [ 0 , 1 ], the random errors ε i are i.i.d. standard normal random variables, and the noise level σ 2 may, or may not, be known. The interest is to recover the function f from the observations y i . Additionally, we assume that the unknown signal f has a bounded L 2 -energy, namely ∫ [ 0 , 1 ] f 2 ( t ) d t, and hence it assumes values from a bounded interval. After applying a linear and orthogonal wavelet transform, the model in (1.2) becomes (see, e.g., [16]), (1.3) c J 0 , k = θ J 0 , k + σ ϵ J 0 , k , k = 0 , … , 2 J 0 − 1 , d j , k = θ j , k + σ ϵ j , k , j = J 0 , … , J − 1 , k = 0 , … , 2 j − 1 , where d j , k ( c j k ), θ j , k and ϵ j , k are the wavelet (scaling) coefficients (at resolution level j and location k) corresponding to y, f and ε, respectively; J 0 and J − 1 are the coarsest and finest levels of detail in the wavelet decomposition. If ϵ’s are i.i.d. standard normal, an arbitrary wavelet coefficient from (1.3) can be modeled as (1.4) [ d | θ ] ∼ N ( θ , σ 2 ) , where, due to (approximate) independence of the coefficients, we omitted the indices j, k. The prior information on the energy bound of the signal energy implies that a wavelet coefficient θ corresponding to the signal part in (1.2) assumes its values in a bounded interval, say Θ = [ − m ( j ) , m ( j ) ], which depends on the level j.

Wavelet shrinkage rules have been extensively studied in the literature, but mostly when no additional information on the parameter space Θ is available. For the implementation of wavelet methods in nonparametric regression problems, we also refer to [4], where some methods are described and numerically compared.

Bayesian shrinkage methods in wavelet domains have received considerable attention in recent years. Depending on a prior, Bayes’ rules are shrinkage rules. The shrinkage process is defined as follows: A shrinkage rule is applied in the wavelet domain and the observed wavelet coefficients d are replaced by with their shrunken versions θ ˆ = δ ( d ). In the subsequent step, by the inverse wavelet transform, coefficients are transformed back to the domain of original data, resulting in a data smoothing. The shape of the particular rule δ ( · ) influences the denoising performance.

Bayesian models on the wavelet coefficients have proved capable of incorporating some prior information about the unknown signal, such as smoothness, periodicity, sparseness, self-similarity and, for some particular bases (e.g., Haar), monotonicity.

The shrinkage is usually achieved by eliciting a single prior distribution π on the space of parameters Θ, and then choosing an estimator θ ˆ = δ ( d ) that minimizes the Bayes risk with respect to the adopted prior.

It is well known that most of the noiseless signals encountered in practical applications have (for each resolution level) empirical distributions of wavelet coefficients centered around zero and peaked at zero. A realistic Bayesian model that takes into account this prior knowledge should consider a prior distribution for which the prior predictive distribution produces a reasonable agreement with observations. A realistic prior distribution on the wavelet coefficient θ is given by (1.5) π = ϵ δ 0 + ( 1 − ϵ ) ξ , where δ 0 is a point mass at zero, ξ is a symmetric distribution on the parameter space Θ and ϵ is a fixed number in [ 0 , 1 ], usually level dependent, that regulates the amount of shrinkage for values of d close to 0. Priors for wavelet coefficients as in (1.5) have been considered by Vidakovic and Ruggeri in [33]. Prior literature on wavelet shrinkage by using Bayesian ideas includes important work by [12, 14, 30, 19], among others. A review by Reményi and Vidakovic in [28] overviews a range of Bayesian wavelet shrinkage methods and provides numerous references on Bayesian wavelet shrinkage strategies.

Our approach here is different from all the prior literature. Here is the motivation.

It is clear that specifying a single prior distribution π on the parameter space Θ can never be done exactly. Indeed the prior knowledge of real phenomena always contains uncertainty and multitude of prior distributions can match the prior belief, meaning that on the basis of the partial knowledge about the signal, it is possible to elicit only a family of plausible priors, Γ. In a robust Bayesian point of view the choice of a particular rule δ should not be influenced by the choice of a particular prior, as long as it is in agreement with our prior belief. Several approaches have been considered for ensuring the robustness of a specific rule, Γ-minimax being one compromise.

In this paper, we incorporate prior information on the boundedness of the energy of the signal (the L 2 -norm of the regression function). The prior information on the energy bound often exists in real-life problems, and it can be modelled by the assumption that the parameter space Θ is bounded. Estimation of a bounded normal mean has been considered in numerous papers, e.g., [20, 7, 11, 15, 9, 17], and in the classic text of Ibragimov and Hasminskii [24].

It is well known that estimating a bounded normal mean represents a challenging problem. In our context, if the structure of the prior (1.5) can be supported by the analysis of the empirical distribution of the wavelet coefficients, even then precise elicitation of the distribution ξ cannot be done without some kind of approximation. Any symmetric distribution supported on the bounded set, say [ − m , m ], can be a possible candidate for ξ.

Let Γ denote the family (1.6) Γ = { π = ϵ δ 0 + ( 1 − ϵ ) ξ , ξ ∈ Γ S [ − m , m ] } , where Γ S [ − m , m ] is the class of all symmetric distributions supported on [ − m , m ], δ 0 is point mass at zero, and ϵ is a fixed constant between 0 and 1. We also require that the distribution ξ does not have an atom at 0.

We consider two models; both assume that wavelet coefficients follow normal distribution (which is a statement about the distribution of the noise), d ∼ N ( θ , σ 2 ). In the Model I, the variance of the noise is assumed to be known, while in the Model II the variance is not known and is given a prior distribution. We have kept treatment of Model II limited. Detailed treatment of Model II would take too much additional space.

The rest of the paper is organized as follows. Section 2 contains mathematical aspects and results concerning the Γ-minimax rules. An exact risk analysis of the rule is discussed in Section 3. Section 4 proposes methods of elicitation of hyper-parameters defining the model. For example, one can imagine “estimating” ϵ and even m itself as hyperparameters by using simulated data. But it is going to depend on the run and none of the theoretical results of our paper would then stand. It can be the topic of an entirely different paper. Performance of the shrinkage rule in the wavelet domain and application to a data set are given in Section 5. In Section 6 we provide the conclusions. Proofs are mostly deferred to Appendix A and B.

In Model II, the variance σ 2 is not known and is given an exponential prior. It is well-known that the exponential distribution is an entropy maximize in the class of all distributions supported on R + with a fixed first moment. This choice is non-informative, in a form of a maxent prior.

The model is: [ d | θ , σ 2 ] ∼ N ( θ , σ 2 ) , [ σ 2 ] ∼ E ( μ ) , f ( σ 2 ) = 1 μ exp − σ 2 μ

The marginal likelihood is double exponential as an exponential scale mixture of normals, [ d | θ ] ∼ DE θ , μ 2 g ( d | θ , μ ) = 1 2 μ exp − 2 μ | d − θ | .

Sketches of proofs of Theorems 1 and 2 are deferred to Appendix A and B.

We compared the performance of the proposed three-point prior estimator with eight existing estimation techniques. The selected existing estimation techniques include: Bayesian adaptive multiresolution shrinker [33] (BAMS), Decompsh [22], block-median and block-mean [1], hybrid version of the block-median procedure [1], blockJS [10], visu-shrink [16], and generalized cross validation [2]. The first five techniques are relying on a Bayesian procedure and they are based on level-dependent shrinkage. The blockJS method uses a level-dependent thresholding, while the visu-shrink and generalized cross-validation techniques use a global thresholding method. Readers can find more details about these techniques in [4].

In the simulation study, we computed the AMSE using the parameter values of l = 6 and k = 2.5 and compared with the AMSE computed for the selected estimation techniques. As can be seen in Fig. 9, the proposed estimator shows comparable and for some signals better performance when compared to the selected estimation methods. In particular, for smooth signals (e.g., wave, heavisine), the three-point prior estimator shows better performance compared to non-smooth signals, such as blip, for instance. Moreover, when comparing the performance of the level-dependent estimation methods, the BAMS estimation method shows competitive (or better) performance for most of the cases. We also investigated the influence of SNR level and the sample size on the performance of proposed estimators, compared to the methods considered. For example, for higher SNR (Fig. 10), the three-point priors-based shrinkage procedure does not provide better performance except for wave, angels, and time shifted sine. In general, the Γ-minimax shrinkage shows comparable or better performance compared to other methods considered, when the SNR is low. Also, for larger sample sizes, the three-point prior shrinkage shows improved performance.