Bayesian Wavelet Analysis Using Non-local Priors

Description

Performs Bayesian wavelet analysis using individual non-local priors as described in Sanyal & Ferreira (2017) <DOI:10.1007/s13571-016-0129-3> and non-local prior mixtures as described in Sanyal (2025) <DOI:10.48550/arXiv.2501.18134>.

Details

The main function is BNLPWA, which has arguments for specifying analysis using individual non-local priors or non-local prior mixtures and various hyperparameter specifications for the wavelet coefficients and scale parameters of the non-local priors. See the manual of BNLPWA for examples.

Author(s)

Nilotpal Sanyal <nsanyal@utep.edu>

Maintainer: Nilotpal Sanyal <nsanyal@utep.edu>

References

Sanyal, Nilotpal. "Nonlocal prior mixture-based Bayesian wavelet regression." arXiv preprint arXiv:2501.18134 (2025).

Sanyal, Nilotpal, and Marco AR Ferreira. "Bayesian wavelet analysis using nonlocal priors with an application to FMRI analysis." Sankhya B 79.2 (2017): 361-388.

Bayesian Non-Local Prior-Based Wavelet Analysis

Description

BNLPWA is the main function of this package that performs Bayesian wavelet analysis using individual non-local priors as described in Sanyal & Ferreira (2017) and non-local prior mixtures as described in Sanyal (2025). It currently works with one-dimensional data. The usage is described below.

Usage

BNLPWA(
  data, 
  func=NULL, 
  method=c("mixture","mom","imom"), 
  mixprob_dist=c("logit","genlogit","hypsec","gennormal"), 
  scale_dist=c("polynom","doubleexp"),
  r=1, 
  nu=1, 
  wave.family="DaubLeAsymm", 
  filter.number=6, 
  bc="periodic"
)

Arguments

Details

Spike-and-slab prior for the wavelet coefficients

For individual MOM prior-based analysis, the spike-and-slab prior for the wavelet coefficient d_{lj} is given by

d_{lj} \mid \gamma_l, \tau_l, \sigma^2, r \sim \gamma_l \; \text{MOM}(\tau_l, \sigma^2, r) + (1-\gamma_l) \; \delta(0),

for individual IMOM prior-based analysis, the spike-and-slab prior for the wavelet coefficient d_{lj} is given by

d_{lj} \mid \gamma_l, \tau_l, \sigma^2, r \sim \gamma_l \; \text{IMOM}(\tau_l, \sigma^2, r) + (1-\gamma_l) \; \delta(0),

and for non-local prior mixture-based analysis, the spike-and-slab prior for the wavelet coefficient d_{lj} is given by

d_{lj} \mid \gamma_l^{(1)}, \gamma_l^{(2)}, \tau_l^{(1)}, \tau_l^{(2)}, \sigma^2, r, \nu \sim \gamma_l^{(1)} \; \text{MOM}(\tau_l^{(1)}, r, \sigma^2) + (1-\gamma_l^{(1)})\gamma_l^{(2)} \;\text{IMOM}(\tau_l^{(2)}, \nu, \sigma^2) + (1-\gamma_l^{(1)})(1-\gamma_l^{(2)}) \; \delta(0),

where the density of the MOM prior is

mom(d_{lj} | \tau_{l}^{(1)},r,\sigma^{2}) = \widetilde{M}_{r} \left(\tau_{l}^{(1)}\sigma^{2}\right)^{-r-1/2} d_{lj}^{2r} \exp\left(-\frac{d_{lj}^{2}}{2\tau_{l}^{(1)}\sigma^{2}}\right), \quad r>1, \tau_{l}^{(1)}>0, \widetilde{M}_{r} = \frac{(2\pi)^{-1/2}}{(2r-1)!!}

and the density of the IMOM prior is

imom(d_{lj} | \tau_{l}^{(2)},\nu,\sigma^{2}) = \frac{\left(\tau_{l}^{(2)}\sigma^{2}\right)^{\nu/2}}{\Gamma(\nu/2)} |d_{lj}|^{-\nu-1} \exp\left( -\frac{\tau_{l}^{(2)}\sigma^{2}}{d_{lj}^{2}} \right),\quad \nu>1, \tau_{l}^{(2)}>0.

Hyperparameter specifications

For non-local prior mixture-based analysis, the available specifications for the mixture probabilities are

1. Logit:

   \gamma_l^{(1)} = \frac{\exp(\theta^{\gamma}_{1} - \theta^{\gamma}_{2}l)} {1 + \exp(\theta^{\gamma}_{1} - \theta^{\gamma}_{2}l)}, \quad \theta^{\gamma}_{1} \in \mathbb{R}, \; \theta^{\gamma}_{2} > 0

   \gamma_l^{(2)} = \frac{\exp(\theta^{\gamma}_{3} - \theta^{\gamma}_{4}l)} {1 + \exp(\theta^{\gamma}_{3} - \theta^{\gamma}_{4}l)}, \quad \theta^{\gamma}_{3} \in \mathbb{R}, \; \theta^{\gamma}_{4} > 0

2. Generalized logit or Richards:

   \gamma_l^{(1)} = \frac{1}{[1 + \exp\{-(\theta^{\gamma}_{1} - \theta^{\gamma}_{2}l)\}]^{\theta^{\gamma}_{3}}}, \quad \theta^{\gamma}_{1} \in \mathbb{R}, \; \theta^{\gamma}_{2},\theta^{\gamma}_{3} > 0

   \gamma_l^{(2)} = \frac{1}{[1 + \exp\{-(\theta^{\gamma}_{4} - \theta^{\gamma}_{5}l)\}]^{\theta^{\gamma}_{6}}}, \quad \theta^{\gamma}_{4} \in \mathbb{R}, \; \theta^{\gamma}_{5},\theta^{\gamma}_{6} > 0;

3. Hyperbolic secant:

   \gamma_l^{(1)} = \frac{2}{\pi} \arctan\left[\exp\left(\frac{\pi}{2} \left(\theta^{\gamma}_{1} - \theta^{\gamma}_{2}l\right)\right)\right], \quad \theta^{\gamma}_{1} \in \mathbb{R}, \; \theta^{\gamma}_{2} > 0

   \gamma_l^{(2)} = \frac{2}{\pi} \arctan\left[\exp\left(\frac{\pi}{2} \left(\theta^{\gamma}_{3} - \theta^{\gamma}_{4}l\right)\right)\right], \quad \theta^{\gamma}_{3} \in \mathbb{R}, \; \theta^{\gamma}_{4} > 0

4. Generalized normal:

   \gamma_l^{(1)} = \frac{1}{2} + \text{sign}(\theta^{\gamma}_{1}-l) \frac{1}{2\Gamma(1/\theta^{\gamma}_{2})} \gamma\left(1/\theta^{\gamma}_{2} ,\left|\frac{\theta^{\gamma}_{1}-l}{\theta^{\gamma}_{3}}\right|^{\theta^{\gamma}_{2}}\right), \quad \theta^{\gamma}_{1} \in \mathbb{R}, \; \theta^{\gamma}_{2},\theta^{\gamma}_{3} > 0

   \gamma_l^{(2)} = \frac{1}{2} + \text{sign}(\theta^{\gamma}_{4}-l) \frac{1}{2\Gamma(1/\theta^{\gamma}_{5})} \gamma\left(1/\theta^{\gamma}_{5} ,\left|\frac{\theta^{\gamma}_{4}-l}{\theta^{\gamma}_{6}}\right|^{\theta^{\gamma}_{5}}\right), \quad \theta^{\gamma}_{4} \in \mathbb{R}, \; \theta^{\gamma}_{5},\theta^{\gamma}_{6} > 0

For individual non-local prior based analysis, gamma_l is defined likewise.

For non-local prior mixture-based analysis, the available specifications for the scale parameters are

1. Polynomial decay:

   \tau_{l}^{(1)} = \theta^{\tau}_{1} l^{-\theta^{\tau}_{2}}, \quad \theta^{\tau}_{1},\theta^{\tau}_{2} > 0

   \tau_{l}^{(2)} = \theta^{\tau}_{3} l^{-\theta^{\tau}_{4}}, \quad \theta^{\tau}_{3},\theta^{\tau}_{4} > 0

2. Double-exponential decay:

   \tau_{l}^{(1)} = \theta^{\tau}_{1} \exp(-\theta^{\tau}_{2} l) + \theta^{\tau}_{3} \exp(-\theta^{\tau}_{4} l), \quad \theta^{\tau}_{1},\theta^{\tau}_{2},\theta^{\tau}_{3},\theta^{\tau}_{4} > 0

   \tau_{l}^{(2)} = \theta^{\tau}_{5} \exp(-\theta^{\tau}_{6} l) + \theta^{\tau}_{7} \exp(-\theta^{\tau}_{8} l), \quad \theta^{\tau}_{5},\theta^{\tau}_{6},\theta^{\tau}_{7},\theta^{\tau}_{8} > 0

For individual non-local prior based analysis, tau_l is defined likewise.

Note: The wavelet computations are performed by using the R package wavethresh.

Value

A list containing the following.

Author(s)

Nilotpal Sanyal <nsanyal@utep.edu>

Maintainer: Nilotpal Sanyal <nsanyal@utep.edu>

References

Sanyal, Nilotpal. "Nonlocal prior mixture-based Bayesian wavelet regression." arXiv preprint arXiv:2501.18134 (2025).

Sanyal, Nilotpal, and Marco AR Ferreira. "Bayesian wavelet analysis using nonlocal priors with an application to FMRI analysis." Sankhya B 79.2 (2017): 361-388.

See Also

wd, wr

Examples

  # Using the well-known Doppler function to 
  # illustrate the use of the function BNLPWA

  # set seed for reproducibility
  set.seed(1)

  # Define the doppler function
  doppler <- function(x) { 
    sqrt(x * (1 - x)) * sin((2 * pi * 1.05) / (x + 0.05)) 
  }

  # Generate true values over a grid of length an integer power of 2 
  n <- 128  # Number of points
  x <- seq(0, 1, length.out = n)
  true_signal <- doppler(x) 

  # Add noise to generate data
  sigma <- 0.2  # Noise level
  y <- true_signal + rnorm(n, mean = 0, sd = sigma)

  # BNLPWA analysis based on MOM prior using logit specification
  # for the mixture probabilities and polynomial decay
  # specification for the scale parameter
  fit_mom <- BNLPWA(data=y, func=true_signal, r=1, wave.family=
    "DaubLeAsymm", filter.number=6, bc="periodic", method="mom", 
    mixprob_dist="logit", scale_dist="polynom")

  plot(y,type="l",col="grey") # plot of data
  lines(fit_mom$func.post.mean,col="blue") # plot of posterior 
  # smoothed estimates
  fit_mom$MSE.mean

  
    # BNLPWA analysis using non-local prior mixtures using generalized 
    # logit (Richard's) specification for the mixture probabilities and 
    # double exponential decay specification for the scale parameter
    fit_mixture <- BNLPWA(data=y, func=true_signal, r=1, nu=1, wave.family=
      "DaubLeAsymm", filter.number=6, bc="periodic", method="mixture", 
      mixprob_dist="genlogit", scale_dist="doubleexp")

    plot(y,type="l",col="grey") # plot of data
    lines(fit_mixture$func.post.mean,col="blue") # plot of posterior 
    # smoothed estimates
    fit_mixture$MSE.mean
  
  
  # Compare with other wavelet methods
  library(wavethresh)
  wd <- wd(y, family="DaubLeAsymm", filter.number=6, bc="periodic")  # Wavelet decomposition  
  
  wd_thresh_universal <- threshold(wd, policy="universal", type="hard")
  fit_universal <- wr(wd_thresh_universal)
  MSE_universal <- mean((true_signal-fit_universal)^2)
  MSE_universal

  wd_thresh_sure <- threshold(wd, policy="sure", type="soft")
  fit_sure <- wr(wd_thresh_sure)
  MSE_sure <- mean((true_signal-fit_sure)^2)
  MSE_sure

  wd_thresh_BayesThresh <- threshold(wd, policy="BayesThresh", type="hard")
  fit_BayesThresh <- wr(wd_thresh_BayesThresh)
  MSE_BayesThresh <- mean((true_signal-fit_BayesThresh)^2)
  MSE_BayesThresh

  wd_thresh_cv <- threshold(wd, policy="cv", type="hard")
  fit_cv <- wr(wd_thresh_cv)
  MSE_cv <- mean((true_signal-fit_cv)^2)
  MSE_cv

  wd_thresh_fdr <- threshold(wd, policy="fdr", type="hard")
  fit_fdr <- wr(wd_thresh_fdr)
  MSE_fdr <- mean((true_signal-fit_fdr)^2)
  MSE_fdr

  # Compare with non-wavelet methods
      # Kernel smoothing
  fit_ksmooth <- ksmooth(x, y, kernel="normal", bandwidth=0.05)
  MSE_ksmooth <- mean((true_signal-fit_ksmooth$y)^2)
  MSE_ksmooth
      # LOESS smoothing
  fit_loess <- loess(y ~ x, span=0.1)  # Adjust span for more or less smoothing
  MSE_loess <- mean((true_signal-predict(fit_loess))^2)
  MSE_loess
      # Cubic spline smoothing
  spline_fit <- smooth.spline(x, y, spar=0.5)  # Adjust spar for smoothness
  MSE_spline <- mean((true_signal-spline_fit$y)^2)
  MSE_spline