## ----include=FALSE------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.align = "center"
)

## ----eval=FALSE---------------------------------------------------------------
# install.packages("wavethresh")
# install.packages("WSwavelet_0.1.0.tar.gz", repos = NULL, type = "source")

## -----------------------------------------------------------------------------
library(WSwavelet)

## -----------------------------------------------------------------------------
set.seed(2026)

n <- 256L
x <- seq(0, 1, length.out = n)

demo_signal <- function(x) {
  1.2 * sin(4 * pi * x) +
    3.0 * exp(-((x - 0.32) / 0.025)^2) +
    2.0 * exp(-((x - 0.73) / 0.05)^2) -
    1.4 * (x > 0.52) +
    1.0 * (x > 0.84)
}

truth <- demo_signal(x)
y <- truth + rnorm(n, sd = 0.55)

plot(
  x, y, type = "l", col = "grey55",
  xlab = "Standardized location", ylab = "Signal value",
  main = "Observed signal and noise-free reference"
)
lines(x, truth, col = "navy", lwd = 2)
legend(
  "topright",
  legend = c("Observed", "Reference"),
  col = c("grey55", "navy"), lty = 1, lwd = c(1, 2),
  bty = "n"
)

## -----------------------------------------------------------------------------
fit_g <- wswavelet(
  y = y,
  likelihood = "gaussian",
  filter.number = 2L,
  family = "DaubExPhase",
  bc = "periodic",
  quadrature_n = 24L
)

plot(
  x, y, type = "l", col = "grey60",
  xlab = "Standardized location", ylab = "Signal value",
  main = "WS-Gaussian reconstruction"
)
lines(x, truth, col = "navy", lwd = 2)
lines(x, fit_g$estimate, col = "firebrick", lwd = 2)
legend(
  "topright",
  legend = c("Observed", "Reference", "WS-Gaussian"),
  col = c("grey60", "navy", "firebrick"),
  lty = 1, lwd = c(1, 2, 2), bty = "n"
)

## -----------------------------------------------------------------------------
fit_g$sigma_hat
fit_g$eta_hat
fit_g$optimization$convergence

## -----------------------------------------------------------------------------
head(
  fit_g$level_summary[, c(
    "level", "coefficients", "spike_probability", "beta", "omega",
    "mean_p0", "mean_pW", "mean_pS"
  )]
)

## -----------------------------------------------------------------------------
names(fit_g$detail[[1L]])

## -----------------------------------------------------------------------------
u <- seq(-1.1, 1.1, length.out = 501L)
plot(
  u, wendland_kernel(u), type = "l", lwd = 2, col = "dodgerblue3",
  ylim = c(0, 1.6), xlab = "Standardized coefficient",
  ylab = "Density", main = "The two standardized slab densities"
)
lines(u, semicircle_kernel(u), lwd = 2, col = "darkorange2")
legend(
  "top", legend = c("Wendland", "Semicircle"),
  col = c("dodgerblue3", "darkorange2"), lty = 1, lwd = 2,
  bty = "n"
)

## -----------------------------------------------------------------------------
curve_g <- evaluate_shrinkage_curve(
  fit_g,
  level_index = 1L,
  d_grid = seq(-3, 3, length.out = 301L)
)

plot(
  curve_g$d, curve_g$estimate, type = "l", lwd = 2,
  col = "firebrick", xlab = "Observed coefficient d",
  ylab = "Posterior mean estimate",
  main = "Fitted coefficientwise shrinkage rule"
)
abline(0, 1, lty = 2, col = "grey40")
abline(h = 0, col = "grey75")
legend(
  "topleft",
  legend = c("Posterior mean", "Identity d"),
  col = c("firebrick", "grey40"), lty = c(1, 2), lwd = c(2, 1),
  bty = "n"
)

## -----------------------------------------------------------------------------
matplot(
  curve_g$d,
  curve_g[, c("posterior_spike", "posterior_wendland",
              "posterior_semicircle")],
  type = "l", lty = 1, lwd = 2,
  col = c("grey25", "dodgerblue3", "darkorange2"),
  xlab = "Observed coefficient d", ylab = "Posterior probability",
  main = "Posterior component probabilities"
)
legend(
  "topright",
  legend = c("Spike", "Wendland", "Semicircle"),
  col = c("grey25", "dodgerblue3", "darkorange2"),
  lty = 1, lwd = 2, bty = "n"
)

## -----------------------------------------------------------------------------
level_fit <- fit_g$detail[[1L]]

post <- level_posterior(
  d = c(-1, 0, 1),
  pi_j = level_fit$spike_probability,
  omega_j = level_fit$omega,
  beta_j = level_fit$beta,
  sigma = fit_g$sigma_hat,
  likelihood = fit_g$likelihood,
  laplace_rate = fit_g$laplace_rate_hat,
  quadrature_n = fit_g$settings$quadrature_n,
  use_exact_wendland = fit_g$settings$use_exact_wendland
)

post$posterior_mean
post$probability

## ----eval=FALSE---------------------------------------------------------------
# fit_l <- wswavelet(
#   y,
#   likelihood = "laplace",
#   filter.number = 2L,
#   family = "DaubExPhase",
#   bc = "periodic",
#   quadrature_n = 24L
# )

## ----eval=FALSE---------------------------------------------------------------
# # One common empirical-Bayes mixture weight:
# fit_constant <- wswavelet(
#   y, likelihood = "gaussian", omega_model = "constant"
# )
# 
# # Wendland-only endpoint:
# fit_wendland <- wswavelet(
#   y, likelihood = "gaussian", fixed_omega = 1
# )
# 
# # Semicircle-only endpoint:
# fit_semicircle <- wswavelet(
#   y, likelihood = "gaussian", fixed_omega = 0
# )

## -----------------------------------------------------------------------------
wavelet_object <- wavethresh::wd(
  y,
  filter.number = 2L,
  family = "DaubExPhase",
  type = "wavelet",
  bc = "periodic",
  verbose = FALSE
)
detail_levels <- 0:(wavethresh::nlevelsWT(wavelet_object) - 1L)

estimate_noise_sd(
  wavelet_object,
  detail_levels = detail_levels,
  mad_levels = 1L
)

## -----------------------------------------------------------------------------
checks <- ws_numerical_checks(quadrature_n = 32L)
head(
  checks[, c(
    "likelihood", "pi", "omega", "max_oddness_error",
    "min_first_difference", "max_support_excess", "all_finite"
  )]
)