A Bayesian wavelet shrinkage rule under LINEX loss function

Abstract

This work proposes a wavelet shrinkage rule under asymmetric LINEX loss function and a mixture of a point mass function at zero and the logistic distribution as prior distribution to the wavelet coefficients in a nonparametric regression model with gaussian error. Underestimation of a significant wavelet coefficient can lead to the bad detection of features of the unknown function, such as peaks, discontinuities, and oscillations. It can also occur under asymmetrically distributed wavelet coefficients. Thus, the proposed rule is suitable when overestimation and underestimation have asymmetric losses. Statistical properties of the rule, such as squared bias, variance, frequentist, and bayesian risks, are obtained. Simulation studies are conducted to evaluate the performance of the rule against standard methods and an application in a real dataset involving infrared spectra is provided.

Keywords:

• Wavelet shrinkage
• LINEX loss
• nonparametric regression
• Bayesian inference
• logistic distribution

1. Introduction

Wavelet shrinkage methods are important tools to be applied in nonparametric regression models under wavelet basis expansion of the unknown function. In the wavelet domain, coefficients of the representation of a function are typically sparse and the few nonzero coefficients are localized at important positions of the function to be recovered, such as peaks, discontinuities, and oscillations. However, the observed (empirical) wavelet coefficients, which are obtained by the application of a discrete wavelet transformation on the data, are not sparse in practice due to the presence of noise. In this sense, wavelet shrinkage rules act by denoising the empirical coefficients and estimating the true wavelet coefficients of the basis expansion of the unknown function. The concept of wavelet shrinkage was introduced in the seminal works of Donoho (1993, 1995), Donoho and Johnstone (1994, 1995), and Donoho et al. (1995).

There are several wavelet shrinkage rules available in the literature, most of them are thresholding. A thresholding rule shrinks a small empirical coefficient that is less than a threshold value to exactly zero. The hard and soft thresholding rules of Donoho and Johnstone (1994) are the most famous and applied in the literature and policies to obtain the threshold value were proposed since then, such as the universal threshold by Donoho and Johnstone (1994), the Stein unbiased risk estimator (SURE) by Donoho and Johnstone (1995), and the cross validation threshold by Nason (1996), among others. Bayesian wavelet shrinkage rules have also been proposed along the last decades. In fact, bayesian procedures allow the incorporation of prior information about the wavelet coefficients in the estimation process such as sparseness level and boundedness (if they are limited) through the elicitation of a prior probabilistic distribution to them. Chipman et al. (1997) proposed a mixture of two zero mean normal distributions as prior, one of them with very small variance, in the same idea of a spike and slab prior. Mixtures of a point mass function at zero and a symmetric and unimodal distribution centered at zero were proposed by Vidakovic and Ruggeri (2001) with the Laplace distribution, the Bickel distribution by Angelini and Vidakovic (2004), the double Weibull distribution by Remenyi and Vidakovic (2015), the logistic distribution by Sousa (2020), and the beta distribution by Sousa et al. (2020).

Although the Bayesian shrinkage rules mentioned above and others are suitable for many real data applications, they were built under symmetric loss functions, in particular, the quadratic loss function, which implies that overestimation and underestimation with the same magnitudes of a given wavelet coefficient have the same loss. Moreover, when the quadratic loss function is considered, mathematical advantages are gained such as differentiation and easier analytical manipulation, in addition to the fact that the Bayesian rule is the posterior expected mean of the wavelet coefficient. However, when overestimation and underestimation should have different losses, the proposed shrinkage rules are not suitable and it is necessary the development of a shrinkage rule that is built under an asymmetric loss function. In this sense, this article proposes a Bayesian shrinkage rule under a point mass function at zero and the logistic distribution as prior to the wavelet coefficient and that is obtained under the so called linear-exponential (LINEX) loss function. This asymmetric loss function that was originally proposed by Varian (1975) penalizes exponentially the underestimation (overestimation) and linearly the overestimation (underestimation) according to convenient choices of its parameters. Under the wavelet estimation point of view, its application is interesting once underestimation of significant wavelet coefficients could not detect relevant features of the unknown function, such as peaks and/or discontinuities for example. Huang (2002) and Torehzadeh and Arashi (2014) proposed policies to obtain the threshold value to the soft thresholding rule under LINEX loss function with a Bayesian interpretation of the estimation process; however, the Bayesian shrinkage rule to be proposed in this work allows prior information about the wavelet coefficients to be taken into account in the Bayesian inference about them.

This article is organized as follows: the statistical model is defined in Section 2. The estimation procedure, including the proposed shrinkage rule under LINEX loss, is described in Section 3. Section 4 is dedicated to simulation studies that were conducted to evaluate the performance of the proposed rule. An application in a real dataset about infrared spectra denoising is done in Section 5. Final considerations are in Section 6.

2. Statistical model

We start with the unidimensional nonparametric regression problem, that is, one observes n = 2^J points $(x_1,y_1),\cdots ,(x_n,y_n)$, J∈ℕ, from the model (2.1) $$\begin{array}{l}{y}_i=f(x_i)+e_i,\end{array}$$(2.1) where $x_1,\cdots ,x_n$ are equidistant scalars, $f\in {\mathbb{L}}^2(\mathbb{R})=\{f:\int f^2<\mathrm{\infty }\}$ is an unknown function and $e_1,\cdots ,e_n$ are independent zero mean normal random errors with unknown common variance σ², σ > 0. The goal is to estimate the function f from the sample $\{(x_i,y_i){\}}_{i=1}^n$ without any assumption regarding the functional structure of f, except that it is squared integrable.

The classical approach to estimate f under the nonparametric point of view is to expand it in terms of some functional basis, such as polynomials, Fourier, splines, and wavelets for example, see Takezawa (2005). In this work, we represent f in model (2.1(2.1) $$\begin{array}{l}{y}_i=f(x_i)+e_i,\end{array}$$(2.1) ) in terms of wavelet basis, (2.2) $$\begin{array}{l}f(x)=\sum _{j,k\in \mathbb{Z}}{\theta }_{j,k}{\psi }_{j,k}(x),\end{array}$$(2.2) where $\{{\psi }_{j,k}(x)=2^{j/2}\psi (2^{j}x-k),j,k\in \mathbb{Z}\}$ is an orthonormal wavelet basis for 𝕃₂(ℝ) constructed by dilations j and translations k of a function ψ called wavelet or mother wavelet and θ_{j, k} are wavelet coefficients that describe features of f at spatial location 2^−jk and scale 2^j or resolution level j. Note that the problem of estimating f becomes a problem of estimating the wavelet coefficients in (2.2(2.2) $$\begin{array}{l}f(x)=\sum _{j,k\in \mathbb{Z}}{\theta }_{j,k}{\psi }_{j,k}(x),\end{array}$$(2.2) ). In fact, wavelet is a function that satisfies the admissibility condition $\int |\mathrm{\Psi }(\omega ){|}^2/|\omega |<\mathrm{\infty }$, where Ψ(⋅) is the Fourier transformation of ψ. The admissibility condition implies that $\int \psi =0$, which motivates the name wavelet. There are several wavelet functions available in the literature, see Vidakovic (1999). For the development of this work, we considered the Daubechies’ compactly supported orthogonal wavelet, that does not have closed form but it is shown in Figure 1.

Figure 1. Daubechies’ compactly supported orthogonal wavelet with eight null moments.

We can rewrite the model (2.1(2.1) $$\begin{array}{l}{y}_i=f(x_i)+e_i,\end{array}$$(2.1) ) in vector notation, (2.3) $$\begin{array}{l}\mathit{y}=\mathit{f}+\mathit{e},\end{array}$$(2.3) where $\mathit{y}=[y_1,\cdots ,y_n{]}^{\prime }$, $\mathit{f}=[f(x_1),\cdots ,f(x_n){]}^{\prime }$, and $\mathit{e}=[e_1,\cdots ,e_n{]}^{\prime }$. We apply a discrete wavelet transform (DWT) on the original data to take them to the wavelet domain. Although DWT is usually performed by fast algorithms such as the pyramidal algorithm, it is possible to represent it by a transformation matrix 𝑾 with dimension n×n which is applied on both sides of (2.3(2.3) $$\begin{array}{l}\mathit{y}=\mathit{f}+\mathit{e},\end{array}$$(2.3) ). Since DWT is linear, we obtain the following model in the wavelet domain (2.4) $$\begin{array}{l}\mathit{d}=\mathit{\theta }+\mathit{\epsilon },\end{array}$$(2.4) where $\mathit{d}=\mathit{W}\mathit{y}=[d_1,\cdots ,d_n{]}^{\prime }$ is the vector of empirical (observed) wavelet coefficients, $\mathit{\theta }=\mathit{W}\mathit{f}=[{\theta }_1,\cdots ,{\theta }_n{]}^{\prime }$ is the sparse vector of unknown wavelet coefficients of f, and $\mathit{\epsilon }=\mathit{W}\mathit{e}=[{\epsilon }_1,\cdots ,{\epsilon }_n{]}^{\prime }$ is the vector of random errors. Thus, we can see the empirical wavelet coefficients 𝒅 as noised versions of the unknown wavelet coefficients 𝜽. Further, the random errors in the wavelet domain ${\epsilon }_1,\cdots ,{\epsilon }_n$ remain independent zero mean normal with common variance σ² since the DWT is orthogonal. We will drop the subindices and consider the model $d=\theta +\epsilon$ for a single wavelet coefficient along the text for simplicity.

The estimation of the vector of wavelet coefficients 𝜽 will be done coefficient by coefficient under a Bayesian framework. In this sense, a prior distribution is assigned to a single wavelet coefficient θ that incorporates our knowledge about its sparsity and symmetry around zero. Then, we consider the prior distribution $\pi (\cdot ;\alpha ,\tau )$ based on a mixture of a point mass function at zero δ₀(⋅) and a symmetric around zero logistic distribution $g(\cdot ;\tau )$ proposed by Sousa (2020), (2.5) $$\begin{array}{l}\pi (\theta ;\alpha ,\tau )=\alpha {\delta }_0(\theta )+(1-\alpha )g(\theta ;\tau ),\end{array}$$(2.5) where $$\begin{array}{l}g(\theta ;\tau )=\frac{\exp \{-\theta /\tau \}}{\tau (1+\exp \{-\theta /\tau \}{)}^2}{\mathbb{I}}_{\mathbb{R}}(\theta ),\end{array}$$

$\alpha \in (0,1)$, τ > 0 and 𝕀_ℝ(⋅) is the indicator function on the real set ℝ. Thus, the prior distribution has two hyperparameters to be elicited, α and τ. As discussed in Sousa (2020), these hyperparameters control the severity of the shrinkage imposed by the Bayesian rule on the empirical coefficients. In fact, the empirical coefficients are shrunk more for higher values of α and smaller values of τ since the prior distribution becomes more concentrated around zero. In the same work, the author suggests values of τ≤10 and the elicitation of α according to Angelini and Vidakovic (2004), (2.6) $$\begin{array}{l}\alpha =\alpha (j)=1-\frac{1}{(j-J_0+1{)}^{\gamma }},\end{array}$$(2.6) where $J_0\le j\le J-1$, J₀ is the primary resolution level, J is the number of resolution levels, J = log₂⁡(n) and γ > 0. They also suggest that in the absence of additional information, γ = 2 can be adopted. Finally, it is necessary to estimate σ for a complete (hyper)parameter elicitation. According to Donoho and Johnstone (1994), based on the fact that much of the noise information present in the data can be obtained on the finer resolution scale, they proposed (2.7) $$\begin{array}{l}\hat{\sigma }=\frac{\text{median}\{|d_{J-1,k}|:k=0,...,2^{J-1}\}}{0.6745}.\end{array}$$(2.7)

3. Estimation procedure

The estimation of the wavelet coefficients vector 𝜽 is done by the application of the Bayesian wavelet shrinkage rule ${\delta }^{\ast }(\cdot )$ associated to the models (2.4(2.4) $$\begin{array}{l}\mathit{d}=\mathit{\theta }+\mathit{\epsilon },\end{array}$$(2.4) ) and (2.5(2.5) $$\begin{array}{l}\pi (\theta ;\alpha ,\tau )=\alpha {\delta }_0(\theta )+(1-\alpha )g(\theta ;\tau ),\end{array}$$(2.5) ) on the empirical wavelet coefficients vector 𝒅, that is, for a single wavelet coefficient θ, its estimator $\hat{\theta }$ is (3.1) $$\begin{array}{l}\hat{\theta }={\delta }^{\ast }(d).\end{array}$$(3.1)

In fact the estimator $\hat{\theta }$ is a shrunk (denoised) version of the empirical coefficient d. Since we have the estimated wavelet coefficients vector $\hat{\mathit{\theta }}=[{\hat{\theta }}_1,\cdots ,{\hat{\theta }}_n{]}^{\prime }$, the function values 𝒇 are then estimated by the inverse discrete wavelet transformation (IDWT) that can be represented by the transpose matrix 𝑾^′, (3.2) $$\begin{array}{l}\hat{\mathit{f}}={\mathit{W}}^{\mathit{\prime }}\hat{\mathit{\theta }}.\end{array}$$(3.2)

The Bayesian rule ${\delta }^{\ast }(\cdot )$ is obtained by minimizing the posterior expected loss ${\mathbb{E}}_{\pi }[L(\delta ,\theta )|d]$, that is (3.3) $$\begin{array}{l}{\delta }^{\ast }(d)=\arg \underset{\delta }{min}{\mathbb{E}}_{\pi }[L(\delta ,\theta )|d],\end{array}$$(3.3) according to a loss function specification L(δ, θ), which is done depending on the estimation problem and its costs. For symmetric losses as example, the well-known quadratic loss function $L(\delta ,\theta )=(\delta -\theta {)}^2$ is generally applied and its associated Bayesian rule is the posterior expected mean ${\delta }^{\ast }(d)={\mathbb{E}}_{\pi }(\theta |d)$, see Robert (2007) for a complete development of Bayesian frameworks. In this work, however, we consider wavelet coefficients estimation problem under asymmetric losses, in the sense that underestimation and overestimation of θ lead to different costs. Thus, the LINEX loss function is a very suitable function to be considered in this context. The next subsections describe the LINEX loss function, the associated Bayesian rule, and its statistical properties.

3.1. LINEX loss function

We consider the unidimensional LINEX loss function by Varian (1975) given by (3.4) $$\begin{array}{l}L(\delta ,\theta )=b[e^{a(\delta -\theta )}-a(\delta -\theta )-1],\end{array}$$(3.4) for a≠0 and b > 0. The sign of a indicates the direction of the exponential and linear losses. In fact, for a > 0, the loss is exponential when δ > θ and linear when δ < θ, that is, the loss is greater in overestimation. On the other hand, for a < 0, the loss is exponential when δ < θ and linear when δ > θ, that is, the loss function penalizes more the underestimation. The parameter b controls the general magnitudes of the losses.

Figures 2(a) and (b) show the LINEX function L(δ, 0) for $a\in \{\pm 0.4,\pm 0.6,\pm 0.8,\pm 1\}$, b = 1 and for a = 1, $b\in \{0.2,0.4,0.6,0.8,1\}$, respectively. Higher absolute values of a imply in higher losses. For example, in Figure 2(a), when a = 0. 8, $L(4,0)=20.33$ while for a = 1, $L(4,0)=49.59$. Thus the loss of an overestimation of 4 units under a = 1 is almost 2.44 times the loss under a = 0. 8. Higher values of b also increase the loss, as can be seen in Figure 2(b). For example, when b = 0. 8, $L(4,0)=39.68$, against 49. 59 when b = 1.

Figure 2. LINEX loss functions L(δ, 0) (i.e., the loss when θ = 0) for $a\in \{\pm 0.4,\pm 0.6,\pm 0.8,\pm 1\}$, b = 1 (a) and for a = 1, $b\in \{0.2,0.4,0.6,0.8,1\}$ (b).

3.2. Bayesian shrinkage rule

The Bayesian rule δ^* is obtained by minimizing the posterior expected loss ${\mathbb{E}}_{\pi }[L(\delta ,\theta )|d]$. According to Zellner (1986), under the LINEX loss function we have (3.5) $$\begin{array}{l}{\mathbb{E}}_{\pi }[L(\delta ,\theta )|d]=b[e^{a\delta }{\mathbb{E}}_{\pi }(e^{-a\theta }|d)-a(\delta -{\mathbb{E}}_{\pi }(\theta |d))-1],\end{array}$$(3.5) and the Bayesian rule is (3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) where log is the natural logarithm. Note that the Bayesian rule under LINEX loss function does not depend on the parameter b. Finally, it is easy to show that under the model, one has that (3.7) $$\begin{array}{ll}& {\mathbb{E}}_{\pi }(e^{-a\theta }|d)\\ & =\frac{\alpha \varphi (d/\sigma )+\sigma (1-\alpha ){\int }_{\mathbb{R}}{e}^{-a(\sigma u+d)}g(\sigma u+d;\tau )\varphi (u)du}{\alpha \varphi (d/\sigma )+\sigma (1-\alpha ){\int }_{\mathbb{R}}g(\sigma u+d;\tau )\varphi (u)du},\end{array}$$(3.7) where ϕ(⋅) is the standard normal density function and the integrals involved in (3.7(3.7) $$\begin{array}{ll}& {\mathbb{E}}_{\pi }(e^{-a\theta }|d)\\ & =\frac{\alpha \varphi (d/\sigma )+\sigma (1-\alpha ){\int }_{\mathbb{R}}{e}^{-a(\sigma u+d)}g(\sigma u+d;\tau )\varphi (u)du}{\alpha \varphi (d/\sigma )+\sigma (1-\alpha ){\int }_{\mathbb{R}}g(\sigma u+d;\tau )\varphi (u)du},\end{array}$$(3.7) ) can be obtained numerically. Figure 3(a) shows the Bayesian rules for σ = 1, α = 0. 9, τ = 1 and $a\in \{\pm 1,\pm 2,\pm 3,\pm 4\}$ on the interval $[-10,10]$. In this figure, it is possible to observe an important feature: the Bayesian rule (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) is not a shrinkage rule in the sense of wavelet shrinkage for all values of a, although the shape of the rules mimics wavelet shrinkage rules. For example, when a = 4, small values of the empirical coefficient d are not shrunk toward zero by the Bayesian rule. Further, if d = 0 then ${\delta }^{\ast }(0)=-0.768$. The main idea of a wavelet shrinkage rule is to shrink small empirical coefficients toward zero, even in an asymmetrically way, as it is the case under asymmetric loss functions or priors for the wavelet coefficients θ.

Figure 3. Bayesian rules (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) for σ = 1, α = 0. 9, τ = 1, and $a\in \{\pm 1,\pm 2,\pm 3,\pm 4\}$ (a) and bayesian shrinkage rules (in wavelet shrinkage sense) for $a\in \{0.5,1,1.5,2,2.5\}$ (b).

Hopefully, we can obtain shrinkage rules for convenient choices of a, as it is shown in Figure 3(b), for $a\in \{0.5,1,1.5,2,2.5\}$. Note that the rules shrink small empirical coefficients toward zero in an asymmetrically manner. Moreover, for empirical coefficients far from zero, the rule acts with different severity according to the sign of the coefficient. Since the LINEX loss penalizes exponentially the overestimation when a > 0, the associated shrinkage rules shrink severely empirical coefficients greater than zero to avoid overestimation of θ, given the prior belief that most of the θ’s are zero, and the severity increases for higher values of a. Although it is not shown in the article, the behavior of the rule is similar when a < 0, but in the opposite sense to avoid underestimation.

3.3. Statistical properties of the shrinkage rule

The squared bias ${\mathrm{Bias}}_{\delta }^2(\theta )=[\mathbb{E}(\delta (d))-\theta {]}^2$ of the shrinkage rules (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) for σ = 1, α = 0. 9, τ = 1 and $a\in \{0.5,1,1.5,2,2.5\}$ were obtained numerically and are shown in Figure 4(a). In fact, the rules have small bias when θ is close to zero and the minimum squared bias is achieved when θ = 0. Since LINEX loss is exponential for δ > θ under positive values of a, that is, the estimation procedure penalizes more overestimation, the squared bias is higher when θ is greater than zero. The shrinkage rule is built under the prior information that most of the wavelet coefficients are zero or close to zero. Then, for coefficients greater than zero, the shrinkage rule tends to underestimate them due to exponential loss imposed by the LINEX loss function against overestimation. Further, the bias increases when a increases, once the loss increases for higher a values. On the other hand, the bias is small when θ is less than zero because the loss is only linear for δ < θ and then, even under the prior belief of sparsity, the estimation procedure allows a good recover of wavelet coefficients less than zero. Analogous interpretation could be made for the frequentist risks $R_{\delta }(\theta )=\mathbb{E}[L(\delta ,\theta )]$ of the rules (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) and displayed in Figure 4(c).

Figure 4. Squared bias (a), variances (b), and frequentist risks (c) of the shrinkage rules (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) for σ = 1, α = 0. 9, τ = 1, and $a\in \{0.5,1,1.5,2,2.5\}$.

The variances ${\mathrm{Var}}_{\delta }(\theta )=\mathbb{E}[(\delta (d)-\mathbb{E}(\delta (d)){)}^2]$ of the shrinkage rules are provided in Figure 4(b). Unlike the squared bias and frequentist risks, the variances have symmetric shapes and are small when θ is close to zero. They increase as $|\theta |$ increases and achieve two maxima. Finally, the variances stabilize close to one for higher values of $|\theta |$. Further, the obtained maxima are higher as the LINEX parameter a decreases.

The general behaviors of the squared bias, variance, and frequentist risk are quite similar. They are small when θ is close to zero, increase until two local maxima and then stabilize for sufficiently high values of $|\theta |$. Moreover, the two local maxima occur at the same locations where their associated shrinkage rules leave the x-axis and start to approximate to the line y = x. For example, when the LINEX parameter a = 1 (red curves in Figure 4), the global maxima of squared bias, variance, and frequentist risk occur at θ = 3. 12 and are equals to 3.77, 1.12, and 4.99, respectively.

Additionally, the squared bias, variance, and frequentist risk of the shrinkage rule (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) seem to achieve limit values for large wavelet coefficients, according to numerical evaluations. Table 1 show those statistical properties for the shrinkage rule with parameters a = 1, σ = 1, τ = 1, and α = 0. 9 applied to some values of θ. When θ increases in the positive sense, ${\mathrm{Bias}}_{\delta }^2(\theta )$ goes to 2.07269, Var_δ(θ) goes to 1.00000, and R_δ(θ) goes to 3.18357. On the other hand, when θ increases in the negative sense, ${\mathrm{Bias}}_{\delta }^2(\theta )$ goes to 0.24752, Var_δ(θ) goes to 1.00000, and R_δ(θ) goes to 1.16713.

Table 1. Squared bias ${\mathrm{Bias}}_{\delta }^2(\theta )$, variances Var_δ(θ), and frequentist risks R_δ(θ) of the shrinkage rule (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) for σ = 1, a = 1, α = 0. 9, and τ = 1.

θ	-100	-50	-10	-5	0	5	10	50	100
${\mathrm{Bias}}_{\delta }^2(\theta )$	0.24752	0.24752	0.24768	0.24294	0.00543	2.47595	2.06677	2.07269	2.07269
Varδ(θ)	1.00000	1.00000	1.00188	1.08033	0.05345	1.43813	0.99589	1.00000	1.00000
Rδ(θ)	1.16713	1.16713	1.16737	1.20755	0.07246	4.05877	3.17286	3.18357	3.18357

Finally, the Bayesian risks r_δ = 𝔼_π[R_δ(θ)] of the shrinkage rules are available in Table 2. We note that the Bayesian risk increases as the LINEX parameter a increases. For example, for a = 0. 5, we obtain r_δ = 0. 1521 and for a = 2. 5 we have r_δ = 0. 4644, which is almost three times the risk for a = 0. 5. Since shrinkage rules under higher values of a shrink in a more severe way empirical coefficients greater than zero (once the loss is exponential), we have that the Bayesian risks associated to more severe shrinkage rules are greater than less severe ones.

Table 2. Bayesian risks of the shrinkage rule (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) for σ = 1, α = 0. 9, τ = 1, and $a\in \{0.5,1,1.5,2,2.5\}$.

a	0.5	1.0	1.5	2.0	2.5
rδ	0.1521	0.2005	0.2478	0.3512	0.4644

3.4. Comparison with quadratic loss function

The quadratic loss function $L(\delta ,\theta )=(\delta -\theta {)}^2$ is perhaps the most considered loss function in statistical decision problems due to several reasons. Under a mathematical point of view, it is more tractably analytically than other loss functions and differentiable. Further, in statistical inference, when it is dealt with an unbiased estimator, the quadratic loss function becomes the variance of the estimator, with several nice properties and a rich theory, see Schervish (1995).

One important feature of the quadratic loss function is the symmetry around θ, which means that underestimation and overestimation with the same magnitudes of a parameter θ by an estimator δ imply in the same polynomial (quadratic) loss. On the other hand, the LINEX loss function penalizes exponentially one side and linearly the other side depending on its parameters, that is, the losses are asymmetric. Figure 5(a) presents the LINEX loss functions for $a=-1$ and a = 1 and the quadratic loss function when θ = 0. Under overestimation, that is, when δ≫0, the LINEX loss for a = 1 is considerably greater than the quadratic loss, as its loss increases exponentially. The LINEX loss for $a=-1$ increases only linearly and has the lowest loss. The same interpretation occurs for the underestimation case, but the roles of the LINEX rules change with each other.

Figure 5. LINEX ($a=-1$, a = 1, and b = 1) and quadratic loss functions for θ = 0, that is, L(δ, 0) (a) and their associated shrinkage rules under the prior model (2.5(2.5) $$\begin{array}{l}\pi (\theta ;\alpha ,\tau )=\alpha {\delta }_0(\theta )+(1-\alpha )g(\theta ;\tau ),\end{array}$$(2.5) ) for α = 0. 9 and $\tau =\sigma =1$ (b).

In the Bayesian framework, the rule under the quadratic loss function is given by the posterior expected value of θ, that is, ${\delta }^{\ast }(d)={\mathbb{E}}_{\pi }(\theta |d)$. In this case, Sousa (2020) showed that under the prior model (2.5(2.5) $$\begin{array}{l}\pi (\theta ;\alpha ,\tau )=\alpha {\delta }_0(\theta )+(1-\alpha )g(\theta ;\tau ),\end{array}$$(2.5) ), the shrinkage rule is given by $$\begin{array}{l}\delta (d)={\mathbb{E}}_{\pi }(\theta |d)=\frac{(1-\alpha ){\int }_{\mathbb{R}}(\sigma u+d)g(\sigma u+d;\tau )\varphi (u)du}{\frac{\alpha }{\sigma }\varphi (\frac{d}{\sigma })+(1-\alpha ){\int }_{\mathbb{R}}g(\sigma u+d;\tau )\varphi (u)du}.\end{array}$$

The shrinkage rules under LINEX loss for $a=-1$ and a = 1 and the shrinkage rule under quadratic loss and prior model (2.5(2.5) $$\begin{array}{l}\pi (\theta ;\alpha ,\tau )=\alpha {\delta }_0(\theta )+(1-\alpha )g(\theta ;\tau ),\end{array}$$(2.5) ) for α = 0. 9 and $\tau =\sigma =1$ are shown in Figure 5(b). In fact, the shrinkage rule under quadratic loss acts symmetrically around zero and shrinks the empirical coefficients with intermediate severity in both sides. The shrinkage rule under LINEX loss with a = 1 shrinks empirical coefficients greater than zero with strong severity and the rule with $a=-1$ is severe for coefficients less than zero.

In general, the choice of the loss function has a great impact on the behavior of the associated shrinkage rule. Symmetric loss functions like the quadratic assigns equals losses for underestimation and overestimation and imply in a symmetric shrinkage rule. Asymmetric loss functions like LINEX are suitable when the loss of underestimation is different to overestimation. In this case, the associated shrinkage rule is then asymmetric around zero.

4. Simulation studies

Simulation studies were conducted to evaluate the performance of the shrinkage rule (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) and to compare it against standard and recent proposed shrinkage and thresholding rules. Wavelet coefficients were generated from a mixture of a point mass function at zero and the beta distribution with right asymmetries given by coefficient of asymmetry (skewness) equals to 0.48 (scenario 1) and 1.73 (scenario 2) and two sample sizes, n = 512 and 2048. Then, empirical wavelet coefficients were obtained from model (2.4(2.4) $$\begin{array}{l}\mathit{d}=\mathit{\theta }+\mathit{\epsilon },\end{array}$$(2.4) ) with iid zero mean gaussian random errors and common variance σ² selected according to two signal to noise ratios (SNR), SNR = 3 and 9. For each replication, r we applied the proposed shrinkage rule under LINEX loss and the comparison rules to the generated empirical coefficients and calculated the mean squared error (MSE)

$$\begin{array}{l}{\mathrm{MSE}}^{(r)}=\frac{1}{n}\sum _{i=1}^n({\hat{\theta }}_i-{\theta }_i{)}^2,\end{array}$$

and the mean absolute error (MAE)

$$\begin{array}{l}{\mathrm{MAE}}^{(r)}=\frac{1}{n}\sum _{i=1}^n|{\hat{\theta }}_i-{\theta }_i|.\end{array}$$

For each scenario of asymmetry, sample size and signal to noise ratio values, 200 replications were run and the averaged mean squared error (AMSE) and the averaged mean absolute error (AMAE) were calculated as performance measures of the rules, given, respectively, by $$\begin{array}{l}\mathrm{AMSE}=\frac{1}{200}\sum _{r=1}^{200}{\mathrm{MSE}}^{(r)},\end{array}$$ and $$\begin{array}{l}\mathrm{AMAE}=\frac{1}{200}\sum _{r=1}^{200}{\mathrm{MAE}}^{(r)}.\end{array}$$

We compared the proposed shrinkage rule (LINEX - Logistic) against the soft thresholding rules with the threshold value chosen according to the universal thresholding (Universal) proposed by Donoho and Johnstone (1994), cross validation (CV) by Nason (1996) and the Stein unbiased risk estimator (SURE) of Donoho and Johnstone (1995), the bayesian shrinkage rules under symmetric beta prior (Beta symmetric) proposed by Sousa et al. (2020) and asymmetric beta prior (Beta asymmetric) by Sousa (2022), and the soft thresholding rule under LINEX loss proposed by Torehzadeh and Arashi (2014). In all the simulated data, the full DWT was applied by using a Daubechies basis with 10 null moments. Further, the adopted primary resolution level was J₀ = 3 for both considered sample sizes. The R package wavethresh (Nason 2022) was used to perform the DWT and to apply the thresholding methods Universal, CV, and SURE.

Table 3 shows the results of the scenario 1 under asymmetric distributed wavelet coefficients with skewness equals to 0.48. The proposed shrinkage rule had great performances for both values of sample sizes and SNR. Actually, the SNR had the main impact on the performances of the rules and the sample size did not have significant influence on them. For SNR = 3, the shrinkage rule under asymmetric beta prior was the best one in terms of AMSE and AMAE. The proposed rule and the shrinkage rule under symmetric beta prior also worked well and had performances close to the asymmetric beta rule one. SURE and the thresholding rule by Torehzadeh and Arashi had intermediate performances and Universal and CV thresholding rules did not worked well. For example, the AMSE and AMAE of the CV thresholding rule were, respectively, 10.24 and 3.38 times the AMSE and AMAE of the asymmetric beta shrinkage rule. For SNR = 9, the proposed rule was the best for both performance measures, with AMSE = 0.0004 and AMAE = 0.0121 when n = 512 and AMSE = 0.0003 and AMAE = 0.0113 when n = 2048. The shrinkage rule under asymmetric beta prior and SURE thresholding rule also had good behavior, but with performances worse than the proposed rule. For example, the SURE thresholding rule (the second best performance) had AMSE and AMAE, respectively, almost 3 and 2 times of the proposed rule ones when n = 2048. On the other hand, the thresholding rule by Torehzadeh and Arashi did not have good performance, which suggests that it works better for low signal to noise ratios.

Table 3. AMSE and AMAE with the standard deviations of the MSEs and MAEs of the shrinkage and thresholding rules for generated wavelet coefficients under asymmetric distribution with skewness equals to 0.48 - scenario 1.

𝒏	Method	𝐒𝐍𝐑 = 𝟑	𝐒𝐍𝐑 = 𝟗
		AMSE(SD) AMAE (SD)	AMSE(SD) AMAE (SD)
512	Universal	0.0304 (0.0018) 0.1141 (0.0043)	0.0038 (0.0002) 0.0416 (0.0013)
	CV	0.0420 (0.0034) 0.1337 (0.0069)	0.0056 (0.0004) 0.0502 (0.0024)
	SURE	0.0090 (0.0009) 0.0652 (0.0041)	0.0010 (0.0001) 0.0224 (0.0014)
	Beta symmetric	0.0045 (0.0005) 0.0428 (0.0036)	0.0049 (0.0008) 0.0338 (0.0003)
	Torehzadeh and Arashi	0.0086 (0.0019) 0.0629 (0.0040)	1.8417 (0.0032) 1.1194 (0.0014)
	LINEX - Logistic	0.0047 (0.0005) 0.0428 (0.0036)	0.0004 (0.0038)0.0121 (0.0010)
	Beta asymmetric	0.0041 (0.0005)0.0396 (0.0036)	0.0007 (0.0014) 0.0245 (0.0093)
2048	Universal	0.0337 (0.0011) 0.1228 (0.0022)	0.0043 (0.0001) 0.0444 (0.0007)
	CV	0.0452 (0.0023) 0.1423 (0.0037)	0.0059 (0.0003) 0.0525 (0.0013)
	SURE	0.0082 (0.0006) 0.0620 (0.0021)	0.0009 (0.0001) 0.0210 (0.0007)
	Beta symmetric	0.0043 (0.0003) 0.0426 (0.0016)	0.0026 (0.0002) 0.0436 (0.0082)
	Torehzadeh and Arashi	0.0076 (0.0005) 0.0622 (0.0018)	1.8517 (0.0018) 1.1231 (0.0047)
	LINEX - Logistic	0.0045 (0.0003) 0.0425 (0.0017)	0.0003 (0.0001)0.0113 (0.0005)
	Beta asymmetric	0.0039 (0.0003)0.0396 (0.0016)	0.0013 (0.0002) 0.0318 (0.0054)

The results of the simulation study under scenario 2 is available in Table 4. In fact, the proposed shrinkage rule had the best performance for both measures in the four combinations of sample sizes and SNR. It suggests that the rule works better for denoising higher asymmetric distributed wavelet coefficients. The symmetric and asymmetric beta shrinkage rules also had good performances and so the SURE thresholding rule. It should be noted that the thresholding rule under LINEX loss by Torehzadeh and Arashi did not work well in these contexts. Further, its performance decreased for SNR = 9 under both sample size values. The standard thresholding rules under Universal and CV policies had intermediate performances relative to the comparison rules, but they were overcome by the proposed rule. For example, under n = 512 and SNR = 3, the AMSE and AMAE of CV rule were about 10.34 and 3.34 times of the proposed rule ones respectively.

Table 4. AMSE and AMAE with the standard deviations of the MSEs and MAEs of the shrinkage and thresholding rules for generated wavelet coefficients under asymmetric distribution with skewness equals to 1.73 - scenario 2.

𝒏	Method	𝐒𝐍𝐑 = 𝟑	𝐒𝐍𝐑 = 𝟗
		AMSE (SD) AMAE (SD)	AMSE (SD) AMAE (SD)
512	Universal	0.1021 (0.0082) 0.2233 (0.0092)	0.0113 (0.0009) 0.0746 (0.0031)
	CV	0.1448 (0.0182) 0.2668 (0.0163)	0.0158 (0.0020) 0.0885 (0.0056)
	SURE	0.0282 (0.0046) 0.1189 (0.0085)	0.0030 (0.0004) 0.0394 (0.0028)
	Beta symmetric	0.0273 (0.0026) 0.1191 (0.0059)	0.0021 (0.0002) 0.0319 (0.0022)
	Torehzadeh and Arashi	0.7657 (0.0074) 0.6411 (0.0087)	0.9632 (0.0057) 0.8099 (0.0074)
	LINEX - Logistic	0.0140 (0.0026)0.0799 (0.0079)	0.0010 (0.0001)0.0208 (0.0020)
	Beta asymmetric	0.0145 (0.0020) 0.0837 (0.0059)	0.0012 (0.0002) 0.0236 (0.0020)
2048	Universal	0.1056 (0.0040) 0.2117 (0.0044)	0.0117 (0.0004) 0.0707 (0.0014)
	CV	0.1451 (0.0086) 0.2492 (0.0076)	0.0160 (0.0009) 0.0828 (0.0026)
	SURE	0.0237 (0.0020) 0.1026 (0.0042)	0.0026 (0.0002) 0.0343 (0.0012)
	Beta symmetric	0.0237 (0.0012) 0.1035 (0.0030)	0.0019 (0.0001) 0.0278 (0.0010)
	Torehzadeh and Arashi	0.7020 (0.0062) 0.5552 (0.0057)	1.0487 (0.0027) 0.8451 (0.0037)
	LINEX - Logistic	0.0115 (0.0012)0.0685 (0.0040)	0.0008 (0.0001)0.0174 (0.0010)
	Beta asymmetric	0.0126 (0.0009) 0.0732 (0.0029)	0.0010 (0.0001) 0.0201 (0.0009)

Finally, Figure 6 presents the boxplots of the MSEs and MAEs of the rules in the scenario 1 (a) and (b) and scenario 2 (c) and (d), for n = 512 and SNR = 3. The boxplots of scenario 1 show the well behavior of the shrinkage rules under beta prior and the proposed rule in both measures. Universal and CV had high MSEs and MAEs relative to the others. On the other hand, the boxplots of scenario 2 emphasize the atypical behavior of the thresholding rule of Torehzadeh and Arashi, with MSEs and MAEs considerably higher than the others ones. Similar behaviors were obtained for the other values of sample size and SNR.

Figure 6. Boxplots of MSEs and MAEs of the shrinkage and thresholding rules of the simulation studies in scenario 1 (a) and (b) and scenario 2 (c) and (d) for n = 512 and SNR = 3. The rules are: 1 - Universal, 2 - CV, 3 - SURE, 4 - Beta symmetric, 5 - Torehzadeh and Arashi, 6 - LINEX - Logistic (proposed rule) and 7 - Beta asymmetric.

Thus, the simulation studies showed that the proposed shrinkage rule can be considered under asymmetric distributed wavelet coefficients context as an alternative to the shrinkage rule under asymmetric beta prior. Their performances were close to each other and for higher kurtosis value, the proposed rule overcome the asymmetric beta shrinkage rule for both measures and all the combinations of sample size and SNR. These results show that although the simulated wavelet coefficients were generated under beta asymmetric distribution with parameters according to the considered kurtosis, which allowed some advantage in the estimation process to the shrinkage rule under asymmetric beta prior in the scenario 1 for SNR = 3, this advantage was not observed in the scenario 2 against the shrinkage rule under LINEX loss and even in the scenario 1 for SNR = 9. In fact, when the chosen parameters of the shrinkage rule under asymmetric beta are not closed enough to the parameters of the distribution which generates the data, its performance is reduced for higher values of kurtosis (scenario 2) and signal-to-noise ratio. The hyperparameters of the asymmetric beta shrinkage rule were chosen according to Sousa (2022) and they are different of the parameters of the beta distribution in the data generation procedure in the simulation studies, which were defined according to the kurtosis value for each scenario.

5. Application in real dataset

Serenoa repens or Saw Palmetto is a palm native from North America whose extracted essential oil is applied in several medical treatments, such as Benign Prostatic Hyperplasia (BPH) in men, disorder of urinary system, sexual impotence, and others. When the essential oil is adulterated by the addition of substances such as olive oil, the contaminated oil loses its chemical and biological properties and becomes inactive against BPH for example. Chemists usually apply infrared (IR) spectroscopy on samples of Saw Palmetto essential oil to detect peaks of absorbance and classify a given sample in pure or adulterated.

Figure 7(a) shows the IR spectra of a sample of pure Saw Palmetto oil in $n=2^{10}=1024$ points. This dataset is available in the R package ChemoSpec by Hanson (2023). Actually the complete dataset consists of IR spectra of 14 samples of pure and adulterated samples in 1868 points. The particular data considered here is the IR spectra of the second oil sample and the considered sample size is the closest dyadic number that is smaller than 1868, i.e, $n=1024=2^{\lfloor {\log }_2(1868)\rfloor }$. Since the goal is to obtain the peaks of absorbance, it is necessary to filter the data and remove noise that can lead to a misclassification of a peak. In this sense, the application of a wavelet shrinkage rule is suitable to denoise this data.

Figure 7. Infrared spectra of Saw Palmetto essential oil (a) and its empirical wavelet coefficients by resolution level after a DWT application (b).

The first step is to obtain the full DWT of the data. Figure 7(b) presents the magnitudes of the empirical wavelet coefficients of the dataset by resolution level. The considered wavelet basis was Daubechies with 10 null moments. In fact large empirical coefficients in the coarser resolution levels are associated to the positions of the absorbance peaks. On the other hand, the empirical coefficients in the finest resolution level are typically related to noise. The wavelet shrinkage rule acts mainly in these coefficients by reducing their magnitudes.

The proposed shrinkage rule (3.6(3.6) $$\begin{array}{l}{\delta }^{\ast }={\delta }^{\ast }(d)=(-1/a)\log ({\mathbb{E}}_{\pi }(e^{-a\theta }|d)),\end{array}$$(3.6) ) was applied in the empirical coefficients and the denoised version of the IR spectra was recover by the IDWT application (3.2(3.2) $$\begin{array}{l}\hat{\mathit{f}}={\mathit{W}}^{\mathit{\prime }}\hat{\mathit{\theta }}.\end{array}$$(3.2) ). The chosen parameter of the LINEX loss function was a = 1, the scale hyperparameter of the logistic density in the prior (2.5(2.5) $$\begin{array}{l}\pi (\theta ;\alpha ,\tau )=\alpha {\delta }_0(\theta )+(1-\alpha )g(\theta ;\tau ),\end{array}$$(2.5) ) was τ = 5, the hyperparameter α was chosen to be level dependent according to (2.6(2.6) $$\begin{array}{l}\alpha =\alpha (j)=1-\frac{1}{(j-J_0+1{)}^{\gamma }},\end{array}$$(2.6) ), the estimated standard deviation was $\hat{\sigma }=0.0004$ according to (2.7(2.7) $$\begin{array}{l}\hat{\sigma }=\frac{\text{median}\{|d_{J-1,k}|:k=0,...,2^{J-1}\}}{0.6745}.\end{array}$$(2.7) ), and the adopted primary resolution level was J₀ = 3. The obtained denoised IR spectra is showed in Figure 8(a). It is possible to observe the action of denoising in the wavenumber interval [500, 700] and [800, 1024], where local noise peaks were smoothed. Figure 8(b) presents a plot of the empirical coefficients in the interval $[-0.10,0.10]$ versus estimated wavelet coefficients. This plot is interesting to see the range of shrunk coefficients, which occurred mainly in the interval of empirical coefficients $[-0.05,0.05]$. Once the signal-to-noise ratio of the data is high, the shrinkage rule was active only in small empirical coefficients typically localized in the finer resolution levels. This procedure can be applied in the other samples for better identification of peaks and noise removal.

Figure 8. Denoised version of the infrared spectra of Saw Palmetto essential oil by the application of the shrinkage rule under LINEX loss function (a) and empirical coefficients in $[-0.10,0.10]$ versus estimated (shrunk) wavelet coefficients (b).

6. Final considerations

The present work proposed a wavelet shrinkage rule under LINEX loss function and a mixture of a point mass function at zero and the logistic distribution as prior distribution to the wavelet coefficients. The application of this rule is suitable when overestimation and underestimation of the coefficients have asymmetric losses. In the wavelet domain, underestimated significant coefficients can lead to a bad detection of important features to be recovered of the unknown function such as peaks, discontinuities and oscillations. In this sense, to assign higher losses to underestimation should be interesting. Further, the simulation studies suggested a good performance of the proposed shrinkage rule in terms of averaged mean squared error and averaged mean absolute error measures to estimate asymmetrically distributed wavelet coefficients against standard shrinkage and thresholding rules.

The magnitudes of the losses and the severity that the empirical coefficients are shrunk by the rule are well controlled according to convenient choices of the parameters of the LINEX function and the hyperparameters of the prior distribution, which is important to have flexibility in real datasets analysis. The impact of the wavelet basis in the estimation process and the proposal of other asymmetric loss functions are suggested as future works.

Disclosure statement

No potential conflict of interest was reported by the author(s).